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INTRODUCTION 

Hyperthermia  and  thermal  ablation  are  thermal  therapies  in  which  the  cytotoxic  effects  of  elevated 
temperatures  in  tissue  are  induced  to  achieve  cell  death  or  render  the  cells  more  vulnerable  to 
ionizing  radiation  and  chemical  toxins.  In  application  of  a  thermal  ablation  approach,  high  intensity 
focused  ultrasound  (HIFU)  and  high  intensity  focused  microwave  (HIFW)  have  been  explored  for 
various  medical  conditions.  Our  research  is  primarily  directed  at  addressing  severe  limitations  of 
HIFU,  that  is,  very  slow  treatment  times  for  large,  deep  lesions,  including  those  close  to  the  chest 
wall,  by  combining  HIFU  and  HIFW,  and  in  each  modality,  introducing  new  techniques  that 
substantially  enhance  the  treatment  results.  Two  modalities  will  capitalize  on  the  advantages  of 
each  in  various  situations,  and  will  allow  comparison  of  those  advantages  as  the  research 
progresses  to  in  vivo  studies.  In  this  collaborative  research,  work  between  the  Microwave  and 
Ultrasound  groups  is  being  conducted  to  develop  a  technique  that  will  synergize  the  fast  heat 
delivery  properties  of  microwaves  with  high  spatial  resolution  of  ultrasound.  This  particular  progress 
report  is  for  the  ultrasound  portion  and  its  future  integration  with  the  microwave  system. 

BODY 

Task  1.  Design  and  construct  an  ultrasound  breast  lesion  treatment  and  targeting  system, 

capable  of  treatment  of  locations  throughout  the  breast  including  next  to  the  chest  wall  and 
compatible  with  microwave  treatment  from  the  other  side. 

Research  Findings:  Progress  on  the  treatment  system  has  included  construction  and  some 
programming  of  the  full  driving  electronics  for  the  298  element  transducer  array  that  is 
compatible  with  treatment  from  both  sides  of  the  compressed  breast  in  the  mammographic 
geometry  and  also  from  outside  a  20  cm  OD  microwave  antenna  array  when  the  array  can  be 
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made  acoustically  transparent  [1],  Fig.  1  shows  the  multichannel  ultrasound  electronics. 
Maximum  acoustic  pressure  achieved  was  approximately  30  MPa  PP.  A  fire  destroyed  almost 
the  entire  laboratory  containing  the  therapeutic  ultrasound  transducer  array  electronics,  but 
insurance  is  paying  for  their  reconstruction  with  some  design  improvements  resulting  from  our 
experience  at  no  additional  cost. 

For  practical  targeting,  we  have  the  GE  Vivid  7  scanner  with  V3  (2000  element)  2D  array 
that  fits  easily  in  the  55  mm  opening  in  the  center  of  the  US  therapy  arrays.  More  advanced 
research  on  ultrasound  imaging  for  targeting  and  correcting  the  microwave  and  US  therapeutic 
and  imaging  beams  was  conducted  in  our  laboratory  with  our  help  by  one  of  Dr.  Moghaddam’s 
students,  Mark  Haynes.  That  work  provides  significant  advances  in  calibration  of  ultrasound 
transducer  arrays  for  full  wave  migration  imaging  that  can  provide  acoustic  properties  of  tissues 
including  speed  of  sound,  scattering,  attenuation  and  compressibility  [2,  3]. 

Task  2.  Evaluate  lesion  production  with  microwave  and  ultrasound  (single  element  and  with 

the  array  transducers)  in  breast  tissue  mimicking  phantoms  at  depths  up  to  5-7  cm. 

Research  Findings:  The  two  modes  are  not  combined  yet.  It  was  determined  that  a 
microwave  antenna  on  a  flat  compression  plate  would  not  give  enough  focal  gain  for  deeper 
heating  in  the  breast  than  at  the  skin.  A  cylindrical  microwave  antenna  was  designed  and 
constructed  as  reported,  with  record  target  to  surrounding  tissue  heating  efficiency.  These 
results  are  described  in  the  report  of  the  initiating  P.l.  While  our  ultrasound  could  be  brought  in 
from  the  bottom  of  the  cylindrical  microwave  antenna,  it  was  deemed  more  effective  to  refine 
the  two  modalities  separately  and  then  initiate  a  new  project  to  allow  ultrasound  transmission 
through  the  microwave  antenna.  Plans  and  quite  satisfactory  progress  with  ultrasound  heating 
and  its  enhancement  have  been  and  are  being  reported  publicly  [1, 4]. 

Task  3.  Test  operation  of  simultaneous  ultrasound  and  microwave  heating  of  tissue-mimicking 

phantoms.  Obtain  initial  data  on  beam  steering  errors  for  validation  of  future  correction  of 
same. 

Research  Findings:  As  discussed  in  Task  2,  simultaneous  operation  has  not  been  achieved 
because  the  ultrasound  transmission  through  the  existing  microwave  substrate  material  is 
unacceptably  high.  However,  experimental  results  using  laboratory  prototypes  generated 
focused  heating  sufficient  for  an  antitumor  effect  deep  in  the  breast  in  tissue-mimicking  gelatin 
phantoms  with  both  microwave  and  ultrasound  separately.  Heating  focal  spot  sizes  were  as 
small  as  1.5  cm  with  microwaves  [5]  and  1.5  mm  with  ultrasound  [1].  Reasonably  acceptable 
lipid/water  emulsions  have  been  found  as  coupling  fluids  for  both  modalities  as  well  as  for 
mammography  [5].  They  currently  leave  an  oily  residue  on  the  skin  [6]. 

Task  4.  Evaluate  future  acceleration  of  acoustic  treatment  by  increasing  sound  absorption 

only  in  the  treatment  zone  by  creation  of  bubbles  in  the  acoustic  focus  using  acoustic  droplet 
vaporization  (ADV);  this  will  also  allow  use  of  higher  microwave  and  US  powers  by  shielding  of 
ultrasonically  distal  skin. 

Research  Findings:  Exceptional  additional  progress  has  been  made  with  acceleration  of  the 
ultrasound  treatment  by  scattering  from  acoustically  vaporized  droplets  [1 , 4],  Emulsions  of 
perfluorocarbon  droplets  (lipid  coated,  C5F12,  0  2.0±0.1  pm,  -99%  <  8  pm  0)  were  used  to 
create  thermal  agents  in  polyacrylamide  phantoms.  The  emulsion  concentration  in  the  gel  was 
3x1 05  droplets  per  mL.  Egg  white  was  incorporated  to  allow  for  visual  inspection  of  the 
phantoms  after  acoustic  exposure.  In  situ  MRI  temperature  monitoring  limited  focal  heating  to 
75 °C.  Lesion  sizes  were  measured  as  a  function  of  applied  acoustic  power.  A  similar,  but 
different  acoustic  array  was  employed  for  this  recent  experiment.  Acoustic  trenches,  as 
sketched  in  Fig.  2,  were  created  by  vaporizing  the  droplets  to  form  bubble  walls  in  a  spiral 
pattern  on  an  underlying  wall  to  trap  sound  scattered  by  the  bubbles  (Fig.  3b).  Exposures  in 
the  middle  of  these  trenches  produced  nearly  uniform  heating  throughout  the  larger  lesion  (Fig. 
3c).  Lesion  volumes  increased  by  a  factor  of  at  least  17  when  comparing  lesion  volumes  in 
phantoms  with  droplets  to  without  droplets.  With  the  use  of  acoustic  droplet  vaporization  (ADV) 
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and  the  resulting  trenches,  a  uniform  ablation  volume  of  15  ml.  was  achieved  in  15  s.  These 
experiments  were  repeated  with  a  commercial  MR  Guided  HIFU  system  with  the  same  results 
[1].  Safety  of  the  emulsions  was  investigated  via  biodistribution  studies  [7]. 

Task  5.  Design  key  features  of  a  clinical  dual  sided  treatment  system  with  degrees  of  freedom 

equivalent  to  those  of  a  mammography  unit  for  treatment  of  the  breast  from  mediolateral, 
through  CC,  to  lateral  medial.  The  design  should  allow  treatment  of  all  breast  tissues 
imageable  mammographically.  The  two-sided  treatment  could  be  either  with  a  single  modality 
or  dual. 

Task  6.  Research  Findings:  Additional  progress  has  been  made  in  planning  two 
approaches  for  the  prone  patient:  1)  a  bowl  microwave  antenna  array  with  the  ultrasound 
transducer  manipulated  around  the  outside  with  a  robotic  arm,  an  initial  design  of  which  was 
shown  in  the  previous  annual  report:  2)  a  cylindrical  microwave  antenna  with  a  circular  ring 
ultrasound  array  for  temperature  measurements  by  speed  of  sound  imaging.  In  this  second 
approach,  after  treatment  of  the  body  of  the  breast,  the  breast  could  be  compressed  against 
the  chest  wall  by  a  membrane  for  chest  wall  heating  by  the  ultrasound  therapeutic  array, 
moved  inside  or  below  the  cylindrical  antenna.  In  both  arrangements,  the  Ultrasound  array 
can  be  split  for  treating  through  the  antenna  bowl  or  cylinder  with  the  flat  side  of  the  array 
close  to  the  chest  wall  when  that  is  required.  These  two  approaches  were  included  in  a  DOD 
Breast  Cancer  Research  Program  Idea  Expansion  proposal  and  the  second  in  an  NIH  ROI 
proposal  on  a  different  body  part.  The  resulting  system  from  the  latter  could  be  used  or 
adapted  for  breast  cancer  treatment.  Neither  proposal  was  funded,  but  resubmissions  are 
planned. 

KEY  RESEARCH  ACCOMPLISHMENTS 

•  Essential  completion  and  calibration  of  an  advanced  298  element  array  system  for  initial 
studies  of  leasioning  for  treatment  from  both  sides  of  the  compressed  breast  or  for  treatment 
through  or  from  below  a  microwave  antenna. 

•  Design  of  methods  to  integrate  ultrasound  and  microwaves  in  a  new,  prone,  uncompressed 
breast  geometry  with  cylindrical  or  bowl  antenna  array. 

•  Further  validation  of  ultrasound  lesion  production  acceleration  with  uniform  coverage  using 
acoustic  droplet  vaporization. 

REPORTABLE  OUTCOMES 

Manuscripts,  publications,  abstracts,  presentations 

All  publications  listed  under  REFERENCES,  below,  except  for  numbers  7  and  8. 

Others: 

Journal  Research  Articles  and  Dissertations  (Copies  of  all  but  the  PhD  dissertation  (Ref.  2)  are 

attached): 

References  2, 3, 5, 6 

M.  Haynes,  J.  Stang,  M.  Moghaddam,  “Microwave  Breast  Imaging  System  Prototype  with  Integrated 

Numerical  Characterization,”  IEEE  International  Journal  of  Biomedical  Imaging,  Volume  2012,  Article 

ID  706365,  1 8  pages,  doi:  1 0.1 1 55/201 2/706365. 

Meeting  Proceedings  Articles  and  Extended  Abstracts: 

Reference  1 . 
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M.  Haynes,  L.  van  Nieuwstadt,  S.  Clarkson,  J.  Stang,  C.  Ward,  M.  Moghaddam,  Ongoing 
Development  of  Microwave  Breast  Imaging  System  Components,  30th  General  Assembly  and 
Scientific  Symposium  of  the  International  Union  of  Radio  Science  (URSI),  Aug  13-20,  2011,  Istanbul 
Turkey  (4  pages). 


Abstracts  and  Meeting  Presentations 
Reference  4. 

J.  Stang,  M.  Haynes,  M.  Moghaddam,  P.  Carson  “Transcutaneous  microwave  thermal  therapy  for 
breast  cancer  treatment  using  image  based  time-reversal  focusing,”  Society  for  Thermal  Medicine 
2012  Annual  Meeting,  April  13-16,  2012,  Portland,  OR. 

J.  Stang,  M.  Haynes,  M.  Moghaddam,  A  Transcutaneous  Microwave  Thermal  Therapy  System 
Prototype  For  Breast  Cancer  Treatment  Using  Image  Based  Time-Reversal  Focusing,  National 
Radio  Science  Meeting,  Jan  4-7,  2012,  Boulder  CO,  invited  abstract. 

M.  Haynes,  J.  Stang,  M.  Moghaddam,  Evaluation  of  a  Full-Cavity  Numerical  Characterization 
Approach  for  an  Experimental  Microwave  Breast  Imaging  System,  National  Radio  Science  Meeting, 
Jan  4-7,  2012,  Boulder  CO,  invited  abstract. 

F.  Padilla,  P.  Carson,  O.D.  Kripfgans,  M.L.  Fabiilli,  A.  Covert,  M.  Zhang,  J.B.  Fowlkes,  J.  Stang,  M. 
AN,  M.  Moghaddam,  "Collaborative  Research  on  Microwave-Ultrasound  Synergistic  Techniques  for 
Treatment  of  Breast  Cancer:  Ultrasound  Ablation  System  and  Method,"  Era  of  Hope  Department  of 
Defense  Breast  Cancer  Research  Program,  Orlando,  Florida  2011). 

M.F.  Zhang  M.  Fabiilli;  P.  Carson,  F.  Padilla,  S.K.  Swanson,  O.  Kripfgans,  J.B.  Fowlkes,  "Spatial 
Control  and  Acceleration  of  Ultrasound  Thermal  Therapy  Using  Acoustic  Droplet  Vaporization," 
Procs.,  Amer.  Inst.  Ultras.  Med.,  J.  Ultras.  Med,  30,  S46,  2011). 

J.  Stang,  P.L.  Carson,  O.D.  Kripfgans,  M.  Fabiilli,  A.  Covert,  M.  Zhang,  J.B.  Fowlkes,  F.  Padilla,  M. 
AN,  M.  Moghaddam,  Ultrasound  Synergistic  Techniques  for  Treatment  of  Breast  Cancer,  201 1  Era  of 
Hope  Conference,  Aug  2-5,  201 1 ,  Orlando  FL. 

J.  Stang,  M.  AN,  O.D.  Kripfgans,  J.B.  Fowlkes,  P.L.  Carson,  M.  Moghaddam,  Microwave  Synergistic 
Techniques  for  Treatment  of  Breast  Cancer,  201 1  Era  of  Hope  Conference,  Aug  2-5,  201 1 ,  Orlando 
FL. 

J.  Stang,  M.  Moghaddam,  O.D.  Kripfgans,  J.B.  Fowlkes,  P.L.  Carson,  Image  Based  Microwave 
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J.  Stang,  M.  Moghaddam,  J.B.  Fowlkes,  P.L.  Carson,  Achieving  Tumoricidal  Dose  with  Focused 
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Patents  and  licenses 

•  University  of  Michigan  Invention  Disclosure  (2011),  File  Number  5273.  Copy  submitted 
separately  by  Initiating  Partner. 

Funding  Applications 

•  NIH  R01  Research  Grant  Application  -  Uniform  Microwave  Hyperthermia  of  Limb  Sarcomas 
with  US  Monitoring  and  Control,  submitted  Feb.  5,  2013. 

•  DoD  CDMRP/BCRP  Idea  Expansion  Award  Application  BC123512P1  -  Integrated  Microwave 
and  Ultrasound  Synergistic  System  for  Treatment  of  Breast  Cancer  12/18/2012. 

Resubmission  of  both  is  planned. 

Employment  and  research  opportunities 

•  Postdoctoral  Research  Fellow,  John  Stang,  received  position  of  Research  Assistant 
Professor  based  on  experience/training  supported  by  this  award. 

•  Mario  Fabiilli,  Postdoctoral  Research  Fellow 

•  Sneha  A.  Nair,  Graduate  Student  Research  Assistant 

•  Oliver  Kripfgans,  Assistant  Research  Scientist 

• 

CONCLUSION 

Others  have  developed  ultrasound  therapy  systems  with  similar  amplitude  and  focusing 
characteristics,  but  not  with  the  split  shell  array  that  can  get  a  large  part  of  the  aperture  close  to  the 
chest  wall  for  treatment  in  the  mammographic  geometry  or  in  the  prone  geometry  in  combination 
with  microwave  treatment.  The  mammographic  geometry  allows  precise,  sure  targeting  by  x-ray  and 
high  resolution  ultrasound  as  verified  by  a  real  time  phased  array  that  can  operate  with  the 
therapeutic  array.  The  prone  geometry  offers  the  promise  of  unprecedented  uniformity  and  other 
control  of  temperature  distributions  by  microwave  and  ultrasound  for  lesion  ablation  as  well  as  whole 
breast  and  adjacent  lymph  nodes  for  local  control.  The  acceleration  of  ultrasound  thermal  and 
ablative  therapy  by  patterned  vaporization  of  IV-injected  perfluorocarbon  microbubbles  is  very 
promising  and  still  unique.  Others  have  begun  working  with  similar  approaches  using  activation  of 
gas  containing  ultrasound  contrast  agents  and  vaporization  of  perfluorocarbon  nanodroplets. 

REFERENCES 

[1 .]  Kripfgans  OD,  Zhang  M,  Fabiilli  M,  Carson  P,  Padilla  F,  Swanson  S,  Mougenout  C,  Fowkes 
B,  "Acceleration  of  ultrasound  thermal  therapy  by  patterned  acoustic  droplet  vaporization,"  J. 
Acoust.  Soc.  Am.  135(1),  537-545,  (2014). 

[2.]  Haynes  MS,  "Full-wave  Nonlinear  Inverse  Scattering  for  Acoustic  and  Electromagnetic 
Breast  Imaging,  Ch.  2,  Acoustic  Inverse  Scattering  Algorithm,  Ch.  8,  Acoustic  Inverse 
Scattering  Experiment,"  PhD  Dissertation,  Univ.  of  Michigan,  201 1 ,  Ann  Arbor,  pp.  281 . 


6 


[3.]  Haynes  M,  Verweij  S,  Moghaddam  M,  Carson  P,  "Transducer  Model  and  Volume  Integral 
Formulation  for  Self-Characterization  of  Commercial  Ultrasound  Probes  in  Transmission 
Acoustic  Inverse  Scattering,"  IEEE  Trans.  Ultras.  Ferroelectrics  Freq.  Control,  accepted  with 
minor  revisions,  (2014). 

[4.]  Fabiilli  ML,  Kripfgans  OD,  Swanson  SD,  Mougenot  C,  Carson  PL,  Zhang  M,  Fowlkes  JB, " 

On  the  Acceleration  of  Ultrasound  Thermal  Therapy  by  Patterned  Acoustic  Droplet 
Vaporization,"  Annual  Meeting  Am.d  Inst.  Ultrasound  Med.,  NY,  NY  (Apr  6-10,  2013). 

[5.]  Stang  J,  Haynes  M,  Carson  P,  Moghaddam  M,  "A  Preclinical  system  prototype  for  focused 
microwave  thermal  therapy  of  the  breast,"  IEEE  Trans  Biomedical  Eng.  59(9),  2431  -  2438, 
(2012). 

[6.]  Carson  PL,  LeCarpentier  GL,  Lashbrook  CR,  Lee  WM,  Goodsitt  MM,  Saitou  K,  Wang  B, 

Baek  S,  Fowkes  B,  Local  Compression  During  Automated  Ultrasound  Scanning  and  Methods 
of  Acoustic  Coupling  (Application  only),  U.S.P.  Office,  No.,  US201 301 1 6570A1  ,May  9,  201 3, 
pp.  25. 

[7.]  Fabiilli  ML,  Piert  MR,  Koeppe  RA,  Sherman  PS,  Quesada  CA,  Kripfgans  OD,  "Assessment  of 
the  biodistribution  of  an  [18F]FDG-loaded  perfluorocarbon  double  emulsion  using  dynamic 
micro-PET  in  rats,"  Contrast  Media  and  Molecular  Imaging  8(4),  366-374,  (2013). 


7 


APPENDICES 


Attach  all  appendices  that  contain  information  that  supplements,  clarifies  or  supports  the  text. 
Examples  include  original  copies  of  journal  articles,  reprints  of  manuscripts  and  abstracts,  a 
curriculum  vitae,  patent  applications,  study  questionnaires,  and  surveys,  etc. 


APPENDIX  1  -  SUPPORTING  DATA 


Figure  1. 

Therapeutic  array  transducer  (black 
arrow)  in  water  tank  aimed  at  needle 
hydrophone  for  automated  output 
calibration  (short  black  arrow).  Amplifier 
and  FPGA  control  boards  and  cables 
(white  arrow)  for  the  298  driver  channels 
are  in  a  clear  plastic  electronics  rack. 


(A)  Backwall  creation 


(B)  Trench  generation 


16  mm 


(C)  Individual  HIFU  exposure 


(D)  Resulting  composite  lesion 


Figure  2.  Diagram  of  an  acoustic 
well  and  first  layer,  spiral  treatment 
pattern. 
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Figure  3.  B-mode  images  of  single  (a)  and  compound  lesions  from  a  US  HIFU  system. 
The  lesion  in  panel  B  is  created  in  a  rapid  succession  of  individual  HIFU  pulses  that 
were  electronically  steered.  Twenty  individual  lesions  were  created  in  500  ms  to  form 
the  spiral  acoustic  trench  (b).  To  only  create  the  bubble  walls,  this  spiral  could  have 
been  at  shorter  duration,  with  thermal  ablation  accomplished  in  the  next  step.  The 
formed  acoustic  trench  was  filled-in  by  20  additional  individual  lesions  which  then 
created  the  solid  lesion  (c). 


APPENDIX  2  -  RESULTING  PUBLICATIONS 

Attached  in  the  following  order: 

[a.]  Kripfgans  OD,  Zhang  M,  Fabiilli  M,  Carson  P,  Padilla  F,  Swanson  S,  Mougenout  C,  Fowkes 
B,  "Acceleration  of  ultrasound  thermal  therapy  by  patterned  acoustic  droplet  vaporization,"  J. 
Acoust.  Soc.  Am.  135(1),  537-545,  (2014). 

[b.]  M.  Haynes,  J.  Stang,  M.  Moghaddam,  “Microwave  Breast  Imaging  System  Prototype  with 
Integrated  Numerical  Characterization,”  IEEE  International  Journal  of  Biomedical  Imaging, 
Volume  201 2,  Article  ID  706365,  1 8  pages,  doi:  1 0.1 1 55/201 2/706365. 

[c.]  Haynes  M,  Verweij  S,  Moghaddam  M,  Carson  P,  "Transducer  Model  and  Volume  Integral 
Formulation  for  Self-Characterization  of  Commercial  Ultrasound  Probes  in  Transmission 
Acoustic  Inverse  Scattering,"  IEEE  Trans.  Ultras.  Ferroelectrics  Freq.  Control,  accepted  with 
minor  revisions,  (2014). 

[d.]  Kripfgans  OD,  Zhang  M,  Fabiilli  M,  Carson  P,  Padilla  F,  Swanson  S,  Mougenout  C,  Fowkes 
B,  "Acceleration  of  ultrasound  thermal  therapy  by  patterned  acoustic  droplet  vaporization,"  J. 
Acoust.  Soc.  Am.  135(1),  537-545,  (2014). 

[e.]  Stang  J,  Haynes  M,  Carson  P,  Moghaddam  M,  "A  Preclinical  system  prototype  for  focused 
microwave  thermal  therapy  of  the  breast,"  IEEE  Trans  Biomedical  Eng.  59(9),  2431  -  2438, 
(201 2).  (Among  the  top  1 8  papers  in  the  journal  in  the  number  of  times  it  was  downloaded  in 
September  2012) 

[f.]  Public  description  of  project  for  BCRP  public  brocure. 


9 


“T 

A  m  i  ] 

V\ 

Ui 

(  ' 

/X  I  |J  Proceedings 

" 

t  / 

|  i  / 

I, H 

IIJ 

7:/,  ■ 

Ultrasound  thermal  ablation  system  and  methods  for  treatment  of  breast  cancer 

O.  D.  Kripfgans,  A.  Covert,  F.  R.  Padilla,  J.  B.  Fowlkes,  M.  L.  Fabiilli,  M.  Moghaddam,  and  P.  L.  Carson 

Citation:  AIP  Conference  Proceedings  1481,  185  (2012);  doi:  10.1063/1.4757332 
View  online:  http://dx.doi.Org/10.1063/1 .4757332 

View  Table  of  Contents:  http://scitation.aip.org/content/aip/proceeding/aipcp/1481  ?ver=pdfcov 
Published  by  the  AIP  Publishing 


Ultrasound  Thermal  Ablation  System  and 
Methods  for  Treatment  of  Breast  Cancer 

O.  D.  Kripfgans,  A.  Covert,  F.  R.  Padilla,  J.  B.  Fowlkes, 

M.  L.  Fabiilli,  M.  Moghaddam3,  P.  L.  Carson 


Departments  of  Radiology  and  aElectrical  Engineering,  University  of  Michigan,  Ann  Arbor,  MI 


Abstract.  The  design  and  initial  testing  of  an  ultrasound  system  capable  of  targeting  and 
treating  breast  cancer  noninvasively  is  presented;  this  system  can  be  used  either  in  the  prone 
position,  thus  compatible  with  MRI  localization  and  microwave-assisted  treatment,  or  in  the 
mammographic  geometry  for  improved  localization.  Two  half-shell  transducers  with  holes  in 
the  centers  for  2D  imaging  arrays  allow  placement  close  to  the  chest  wall  around,  or  on  either 
side  of,  the  breast.  To  drive  the  pair  of  149  element  arrays,  we  are  building  and  programming 
circuitry  to  output  radiofrequency  (RF)  tone  bursts  at  1.5  MHz.  The  low  frequency  may  allow 
aperture  sharing  with  microwave  antennas  and  yet  is  high  enough  for  acoustic  droplet 
vaporization  at  4  MPa  peak  rarefactional  pressure  for  pulse  trains  of  2  s  at  1%  duty  factor  and  far 
less  for  CW  exposures.  Five  Altera  DEI  development  boards  individually  output  RF  tone  bursts 
through  64  programmable  channels  each.  A  MATLAB -based,  graphical  user  interface  controls 
the  system.  Acoustic  output  tests  at  low  pressures  on  a  few  array  elements  and  electrical 
characterization  of  all  elements  suggest  the  necessary  pressures  can  be  achieved  with  either  half 
of  the  full  0.7  f-number  array  using  inexpensive  amplifiers. 

Keywords:  HIFU,  Minimally  Invasive  Therapy,  Transducer  Array,  FPGA 
PACS:  43.80.Vj,  43.80.Sh 


INTRODUCTION 

This  research  is  directed  primarily  toward  eliminating  the  most  severe  limitation  of 
high  intensity  focused  ultrasound  (HIFU),  that  is,  very  slow,  and  occasionally 
incomplete  treatments,  particularly  for  large,  deep  lesions,  including  those  close  to  the 
chest  wall.  When  inhomogeneous  ultrasound  attenuation  or  small,  diffuse  metastases 
are  expected,  combining  HIFU  and  high  intensity  focused  microwave  (HIFW)  may  be 
advantageous.  The  use  of  two  modalities  might  not  be  much  more  complex  clinically 
than  treatment  with  one  modality.  Here  we  address  the  ultrasound  system  and  its 
potential  operation  in  the  mammographic  geometry.  We  are  gaining  confidence  that 
whole  breast  ultrasound  imaging  can  be  performed  in  this  geometry  with  the  images 
closely  aligned  to  those  from  3D  mammography  (digital  breast  tomosynthesis)  [1,  2], 
Coupling  to  the  breast  at  a  much  lower  frequency  for  HIFU  and  imaging  the  treatment 
with  a  2D  array  should  allow  placement  of  the  focus  anywhere  in  the  breast  [3], 
including  over  tumors  close  to  the  chest  wall,  without  substantial  chest  wall  heating 
Figure  1).  We  are  also  making  progress  with  several  other  imaging  modalities  that  are 
compatible  with  this  geometry  and  capable  of  improving  guidance  and  assessment  of 
treatment,  including  microwave,  speed  of  sound  [4]  and  vascularity  imaging  [5,6], 
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The  development  and  testing  of  HIFU  phased  arrays  is  of  great  interest  due  to  their 
therapeutic  potential.  While  tissue  ablation  by  HIFU  has  already  gained  approval  by 
the  FDA  for  the  treatment  of  uterine  fibroids,  development  continues  on  other  possible 
applications  that  are  less  forgiving  of  incomplete  treatment,  such  as  thermal  necrosis 
of  malignant  masses.  Researchers  are  investigating  how  to  maintain  the  focus  of  an 
array  through  or  around  aberrating  structures,  minimize  damage  to  tissue  surrounding 
the  target  area,  safely  increase  focal  power,  and  steer  the  focus  in  three-dimensional 
space. 

Since  HIFU  can  heat  a  relatively  small  volume  of  tissue  at  a  time,  this  leads  to  long 
treatment  times  that  are  limited  by  the  heating  of  overlying  tissues.  Ongoing  work 
with  microbubbles  and  HIFU  has  demonstrated  the  ability  to  greatly  increase  the  local 
absorption  of  sound  in  the  treated  tissue  volume  with  less  absorption  of  overlying 
tissues  in  the  beam  paths.  These  microbubbles  can  be  1)  created  by  very  high 
intensity  ultrasound;  2)  injected  as  ultrasound  contrast  agents;  3)  generated  by  acoustic 
vaporization,  as  planned  in  this  work,  of  injected  perfluorocarbon  droplets  only  in  the 
transducer  focal  zone  [7], 

In  order  to  perform  these  techniques  with  a  phased  array,  each  element  must  have 
its  own  fully  programmable  driver  channel.  This  allows  the  output  from  each  element 
to  be  phased  and  adjusted.  Each  channel  must  have  a  source  that  generates  the 
waveform,  an  amplification  device,  and  in  most  cases,  a  matching  circuit  to  raise  the 
impedance  of  the  driver  circuit  to  equal  that  of  the  element. 


2D  imaging 
array 


FIGURE  1.  Schematic  of  a  breast  in  the  mammographic  position  with  a  therapeutic  array  above  and 
below  tissue  to  be  treated.  The  use  of  half-shells  allows  for  focusing  and  treatment  close  to  the  chest 

wall. 
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Similar  projects  utilize  such  systems  to  drive  their  arrays.  Dr.  Kullervo  Hynynen's 
1372-element,  concentric  array  is  driven  by  an  in-house  developed  RF  system  with 
2000  individual  channels  [8,  9].  This  system  is  controlled  by  computer  via  a  high 
speed  digital  interface  (National  Instruments'  PCI-6534).  Dr.  Mathias  Fink's  quasi¬ 
random,  300-element  array  is  driven  by  a  300-channel  electronic  driving  system  [10]. 
Three-hundred  matching  circuits  are  used  to  raise  the  impedance  of  the  circuits  to  the 
impedance  of  the  elements. 

This  paper  describes  the  design,  programming,  fabrication,  and  evaluation  of  a  RF 
driver  system  for  a  298-element,  concentric  array  (Imasonic  SA,  Besancon,  France) 
composed  of  inexpensive  Field-Programmable  Gate  Arrays  (FPGA)  development 
boards,  amplifier  boards,  and  matching  circuitry.  The  system  was  designed  with  the 
following  characteristics: 

1 .  Fully  programmable  individual  broadband  channels  capable  of  generating  tone 
burst  digital  waveforms  with  10  ns  resolution 

2.  Frequency  range  of  1.5  kHz  to  50  MHz 

3.  Up  to  20  W  per  channel 

4.  Conducive  Graphical  User  Interface  (GUI) 

MATERIALS  AND  METHODS 
Therapeutic  Array  Design 

A  high-power,  split-array  design  was  pursued  for  combination  of  microwave  and 
ultrasonic  heating.  Active  elements  range  from  55  mm  diameter  (inner  ring)  to  130 
mm  diameter  (outer  ring).  Figure  2  illustrates  the  split  between  the  two  half-shells 
resulting  in  149  aperture  elements  each.  Separability  facilitates  the  incorporation  of  an 
additional  microwave-based  treatment  device.  Mechanical  focusing  is  achieved  by  a 
spherical  radius  of  curvature  of  100  mm,  yielding  an  f-number  of  0.76.  The  elements 
are  divided  over  seven  rings  spaced  approximately  5.7  mm  apart  in  a  concentric 
fashion. 


Signal  Generation 

Digital  waveforms  are  generated  by  off-the-shelf  FPGA  control  boards. 
Specifically,  6  development  boards  (DEI,  Altera  Corporation,  San  Jose,  CA,  Figure  3) 
were  employed  for  triggering,  synchronization,  and  signal  generation.  Software  was 
developed  in  Yerilog  by  use  of  Quartus  II. 

Each  FPGA  board  can  control  a  maximum  of  64  ultrasound  channels.  Each 
channel  can  be  controlled  remotely  for  center  frequency,  number  of  transmitted 
waveform  cycles,  waveform  delay  for  additional  electronic  focusing,  and  duty-factor 
for  pulse-width-modulation  waveform  amplitude  shaping.  These  values  are 
communicated  to  the  board  via  read/write  USB  (universal  serial  bus)  from  a  control 
computer.  Transmission  is  triggered  either  by  software  control  or  an  external  TTL 
input  (3.3  Vpp,  >50  ns,  positive  edge).  For  use  in  external  devices  such  as  an 
oscilloscope,  a  trigger-out  signal  is  provided  (3.3  Vpp,  2  ps,  positive  edge). 
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Two  additional  user  programmable  channels  are  available.  A  master/client  setup 
synchronizes  both  triggering  as  well  as  FPGA  clock  cycles.  One  master  board  is 
employed  to  provide  a  common  100  MHz  digital  clock  signal  for  the  signal  generating 
client  boards  as  well  as  for  triggering.  Waveform-delay  as  well  as  transmit  frequency 
selection  is  bound  to  10  ns  resolution.  Sixteen  bit  counters  allow  for  up  to  650  ps 
time-delay  between  channels  as  well  as  65,000  cycle  waveforms.  Very  long  HIFU 
pulse  trains  are  accommodated  by  repeated  trigger  events.  This  allows  for  up  to  85  s 
quasi  continuous  excitation. 


-60  —40  -20  0  20  40  60 

sc  [irm] 

(a)  (b) 

FIGURE  2.  (a)  Design  photos  of  the  therapeutic  array  transducer.  Each  half-shell 
contains  7  concentric  rings  with  a  total  of  149  transducer  elements,  (b)  A  broad  focus 
at  a  distance  of  10  cm  is  achieved  by  shaping  the  array  in  an  elliptical  fashion  as 
shown  here  (Note:  each  element  consists  of  two  triangular  picture  elements.) 
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FIGURE  3.  Employed  off-the-shelf  FPGA  development  hardware.  PC  to  FPGA 
communication  is  facilitated  by  the  use  of  an  USB  protocol.  Digital  waveforms  are 
transmitted  to  the  amplifier  electronics  via  the  general  purpose  I/O  expansion  headers. 
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User  Interface 


Two-way  USB  communication  between  a  host  computer  and  all  6  FPGA  boards  is 
facilitated  by  an  in-house  DLL  driver  and  MATLAB  (MathWorks,  Natick  MA).  A 
conducive  graphical  user  interface  allows  for  channel-by-channel  waveform  control  as 
described  above,  triggering,  and  motorized  3D  positioning  of  the  array  transducer 
(6K4  Motion  Controller,  Parker  Hannifin  Corporation,  Cleveland,  Ohio). 

Acoustic  Measurements 

The  acoustic  output  was  sensed  with  a  needle  hydrophone  (HNC-1500,  Onda 
Corporation,  Sunnyvale,  CA)  and  captured  using  a  digital  oscilloscope  (WaveJet 
324A,  LeCroy  Corporation,  Chestnut  Ridge,  NY). 

RESULTS 

Non-amplified  signal  generation  produced  a  60  dB  signal-to-noise  waveform  with  a 
200  kPa  peak-to-peak  pressure,  for  a  3.3  Ypp  input  of  a  10  cycle  tone  burst  (1.5  MHz 
center  frequency).  The  generated  RF  waveforms  show  a  root-mean-square  jitter  of  2.5 
ns  with  a  standard  deviation  of  1.9  ns,  well  under  the  inherent  10  ns  system  clock 
resolution. 


Peak-Output  Performance 

Short-time  performance  allowed  for  up  to  18  Watts  of  electrical  power  per  channel 
and  a  usual  efficiency  of  70%.  Typically,  15-cycle  tone  bursts  were  transmitted  with 
up  to  1  kHz  pulse  repetition  rates.  The  single-channel  performance  of  amplified  rf- 
signals  yielded  130  kPa  peak-to-peak  in  the  natural  focus  of  the  array.  Assuming 
linear  superposition  of  298  elements,  this  would  correspond  to  approximately  40  MPa 
peak-to-peak,  not  taking  nonlinear  propagation  into  account. 

Time- Average  Power  Performance 

Long-time  performance  was  tested  for  up  to  1  Watt  of  electrical  power  per  channel. 
Trigger  signals  were  employed  to  ensure  >98%  duty  cycle  tone-burst  output  for  5  s. 
Single-channel  performance  yielded  17.5  kPa  peak-to-peak  acoustic  output.  Again, 
assuming  linear  superposition  of  298  elements,  this  would  correspond  to  >5  MPa 
peak-to-peak  pressure  wave. 


CONCLUSIONS 

In  general,  HIFU  arrays  are  not  innovative,  though  a  split  ring  design  is  new  and 
will  be  useful  for  operation  close  to  the  chest  wall.  The  driving  circuits  employ 
standard,  non-expensive,  off-the-shelf  components.  Commercial  signal  generators  that 
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would  meet  the  aforementioned  requirements  would  likely  be  more  expensive. 
National  Instrument’s  PCI-6542  digital  waveform  generator  provides  32  digital  output 
channels  with  an  inter-channel  jitter  of  0.6  ns  for  $4999  [9].  The  FPGA  design 
described  in  this  paper  will  be  1/10  of  that  cost.  The  obtained  acoustic  output 
pressures  are  acceptable  for  both  HIFU  as  well  as  hyperthermia  mode  excitation. 
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The  increasing  number  of  experimental  microwave  breast  imaging  systems  and  the  need  to  properly  model  them  have  motivated 
our  development  of  an  integrated  numerical  characterization  technique.  We  use  Ansoft  HFSS  and  a  formalism  we  developed 
previously  to  numerically  characterize  an  S-parameter-  based  breast  imaging  system  and  link  it  to  an  inverse  scattering  algorithm. 
We  show  successful  reconstructions  of  simple  test  objects  using  synthetic  and  experimental  data.  We  demonstrate  the  sensitivity 
of  image  reconstructions  to  knowledge  of  the  background  dielectric  properties  and  show  the  limits  of  the  current  model. 


1.  Introduction 

A  number  of  experimental  systems  for  microwave  breast 
imaging  have  been  developed  in  recent  years.  These  systems 
test  full-wave  inverse  scattering  algorithms  [1-4]  as  well 
as  synthetic  aperture  beam  focusing  techniques  [5].  While 
imaging  algorithms  abound  in  the  literature,  techniques  to 
properly  model,  characterize,  and  calibrate  these  systems 
have  lagged  behind  algorithm  development.  Investigators 
have  started  to  identify  characterization  as  a  major  task, 
which  must  be  addressed  in  order  to  fully  evaluate  the  effi¬ 
cacy  of  microwave  imaging  for  breast  cancer  detection.  Part 
of  this  evaluation  involves  separating  modeling  errors  from 
intrinsic  algorithm  artifacts  in  the  final  images.  Thus,  there  is 
a  need  for  accurate  models  of  experimental  systems,  as  well 
as  methods  that  efficiently  incorporate  these  models  into  the 
imaging  algorithms. 

The  task  of  characterizing  a  microwave  breast  imaging 
system  for  inverse  scattering,  as  compared  to  a  free- space 
system,  is  complicated  by  several  factors.  Specifically,  the 
antennas  are  not  isolated  in  the  background  media  but 
exist  as  part  of  the  surrounding  structure.  Also,  compact 
arrangements  of  many  antennas  create  a  cavity-like  imaging 
geometry,  and  the  transmitter  incident  fields  include  all 
background  multiple  scattering.  Finally,  the  antennas  and 


object  are  in  each  others  near- fields,  so  object- cavity  scatter¬ 
ing  should  be  modeled. 

In  trying  to  characterize  breast  imaging  systems,  investi¬ 
gators  have  turned  to  full  numerical  simulation.  The  antenna 
cavity  in  [6]  was  modeled  using  Ansoft  HFSS  and  only 
used  for  antenna  design  and  sensitivity  analysis.  In  [7], 
dipole  sources  of  an  inverse  scattering  experiment  were 
modeled  with  HFSS  and  calibration  constants  used  to  scale 
the  antenna  incident  fields.  HFSS  has  also  been  used  to 
obtain  antenna  incident  fields  in  a  near- field  and  open, 
antenna  setup  [8];  however,  ad  hoc  methods  have  been 
used  to  calibrate  the  scattered  field  S-parameter  data  for 
the  inverse  scattering  algorithm.  In  more  recent  work  [9], 
CST  Microwave  Studio  was  used  to  study  and  tune  antenna 
performance  in  a  breast  imaging  cavity.  Also,  finite-volume 
time-domain  solvers  of  [10]  modeled  wide-band  antennas 
for  time-domain  beam  focusing.  The  most  complete  work 
to  date  is  [11],  where  an  FEM  forward  solver  is  used  to 
simulate  the  entire  breast  in  the  presence  of  the  antennas, 
but  computational  complexity  remains  a  challenge.  Despite 
the  growing  use  of  numerical  solvers  to  model  breast 
imaging  systems,  there  has  been  no  clear  or  formal  way  of 
incorporating  the  results  from  full-wave  numerical  models 
into  the  imaging  algorithms. 
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Figure  1:  Microwave  network  model  of  cavity  and  scattering  object. 
S-parameters  are  measured  between  the  reference  planes  on  the 
transmission  lines. 


Figure  2:  Breast  imaging  system  prototype.  The  imaging  cavity  is 
connected  to  the  VNA  through  a  solid-state  switching  matrix.  A 
rotator  is  mounted  above  and  turns  suspended  objects  for  multiple 
transmitter  views. 


The  task  of  characterizing  any  inverse  scattering  system 
can  be  divided  into  three  parts:  (1)  determining  the  incident 
fields  produced  by  the  antennas  in  the  absence  of  the  object, 
(2)  determining  the  background  dyadic  Greens  function, 
that  is,  modeling  the  interactions  between  the  object  and 
its  surroundings  if  necessary,  and  (3)  linking  the  volume 
integrals  in  the  imaging  algorithms  to  measurable  transmit 
and  receive  voltages.  The  purpose  of  this  paper  is  to  show 
how  we  use  FFFSS  and  a  formalism  we  developed  in  previous 
work  [12]  to  solve  parts  (1)  and  (3)  of  this  characterization 
problem,  in  order  to  make  a  numerical  characterization  and 
inverse  scattering  algorithm  consistent  with  an  5-parameter 
based  prototype  breast  imaging  system. 

The  inverse  scattering  algorithm  we  use  is  the  Born  it¬ 
erative  method  (BIM)  with  multivariate -covariance  cost 
function  [13-15].  This  cost  function  allows  us  to  exper¬ 
imentally  choose  the  regularization  parameters  based  on 
our  prior  knowledge  of  system  noise  and  expected  range  of 
permittivities.  The  forward  solver  used  in  the  BIM  requires 


Figure  3:  Imaging  cavity.  Twelve  panels  with  three  bow-tie  antennas 
each  are  solder  together  and  to  a  conducting  plate. 


Figure  4:  HFSS  CAD  model  of  the  imaging  cavity.  Twelve  panels 
contain  three  bow-tie  antennas  each.  The  bottom  of  the  cavity  is 
PEC,  it  is  filled  with  the  coupling  fluid  up  to  the  visible  line,  and  the 
top  surface  radiates  to  air. 


the  background  dyadic  Greens  function  and  finding  it  con¬ 
stitutes  part  (2)  of  the  characterization  problem  mentioned 
above.  For  convenience  we  use  the  lossy  free-space  dyadic 
Greens  function  and  give  some  numerical  and  experimental 
justification  for  this.  Fully  modeling  the  multiple  scattering 
between  the  breast  and  the  imaging  structure  in  the  forward 
solver  is  not  trivial  and  we  discuss  it  in  The  Appendix. 

We  validate  our  methods  with  a  combination  of  sim¬ 
ulation  and  experiment.  We  first  present  the  formalism  of 
[12]  in  the  context  of  cavity  problems.  We  then  explain  our 
experimental  setup,  which  consists  of  a  cylindrical  imaging 
cavity  with  printed  antennas,  solid-state  switching  matrix, 
and  water/oil  coupling  medium.  The  FFFSS  numerical  model 
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Figure  5:  Measured  and  simulated  magnitude  and  phase  of  incident  S2i  between  each  of  the  three  transmitting  antennas  and  all  receivers. 
Solid:  measured.  Dots:  HFSS.  The  groupings  from  left  to  right  are  the  eleven  receivers  of  each  level  (middle,  top,  and  bottom),  repeated  for 
the  three  transmitters  (middle,  top,  and  bottom),  plotted  counterclockwise  when  viewed  from  above  for  a  given  receiver  level.  For  example, 
data  38:48  are  middle  receivers  and  top  transmitter.  The  magnitude  and  phase  agree  best  for  transmitters  and  receivers  on  the  same  level 
(i.e.,  data  1 : 1 1,  50 : 60,  and  98  : 108). 


0  100  200 
(mm) 

Figure  6:  HFSS  CAD  model  of  the  imaging  cavity  with  mesh  of 
unassigned  sheets  to  constrain  the  adapting  meshing  of  HFSS  for 
field  interpolation.  Sheets  are  spaced  every  5  mm  in  each  direction. 


is  presented  and  the  simulation  results  are  compared  to 
those  of  experiment.  We  form  3D  images  of  the  relative 
permittivity  and  conductivity  using  both  HFSS  synthetic 
data  and  experimental  data  for  simple  targets.  We  also 
present  findings  on  the  sensitivity  of  image  reconstructions 
to  the  accuracy  of  modeling  the  background  electrical 
properties. 

Future  work  includes  continuing  the  validation  of  our 
methodology,  experimentally  imaging  more  realistic  breast 
phantoms,  designing  a  hemispherical  imaging  cavity,  inves¬ 
tigating  practical  solutions  to  modeling  the  breast- cavity 
scattering  interactions,  and  developing  a  clinical  imaging 
system. 
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Figure  7:  HFSS  convergence  with  number  of  tetrahedra  for  each 
adaptive  meshing  step. 


2.  Formulation  with  Source  Characterization 

2.1.  Traditional  Volume  Integral  Equations.  The  electric  field 
volume  integral  equation  (VIE)  for  an  inhomogeneous 
distribution  of  permittivity  and  conductivity  is  given  by 

E(r)  =  Einc(r)  +  A^  jG(r,r')  ■  (se(r')  +  ^E(r  )dV\ 

(1) 

where  E(r)  and  Einc(r)  are  the  total  and  incident  fields, 
respectively,  and  r  is  the  position  vector.  The  lossless 
background  wave  number  is  given  by  k* 1 20  =  co2q0eb,  where 
the  background  permittivity  is  €b  =  £0£rb  with  relative 
permittivity  erb.  The  object  contrast  functions  are  defined: 

eh8e(r)  =  e(r)  -  eh, 

(2) 

d<r(r)  =  o{ r)  -  oh , 

where  Ob  is  the  background  conductivity.  The  quantity 
Se( r)  is  unitless  and  5cr(r)  is  an  absolute  measure  of 
conductivity  with  units  of  Siemens  per  meter  and  G(r,r')  is 
the  background  dyadic  Greens  function. 
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Figure  8:  Crosscuts  through  the  center  of  the  cavity  of  the  z-component  of  the  incident  electric  field  due  to  the  middle  transmitter.  The  scale 
is  20logl0(  |  Re{£Z)inC}  I )  of  the  unnormalized  field,  (a)  Horizontal  x-y,  and  (b)  vertical  x-z,  and  (c)  vertical  y-z  planes. 


Defining  the  scattered  field  as 

Esca(r)  =  E(r)  -  Einc(r)  (3) 

and  restricting  the  observation  point  r  to  points  outside  the 
object  region  in  (1),  we  can  write  the  VIE  for  the  scattered 
field  concisely  as 

Esca(r)  =  |  G(r,  r')  ■  0(r')E(r')dV' ,  (4) 

where  we  define  the  following  object  function: 

0(r)  =  k20(se(r)  +  i^^-\  (5) 

V  €bC0  ) 

In  the  context  of  inverse  scattering,  (1)  represents  the 
solution  to  the  wave  equation  in  the  object  domain,  while  (4) 
relates  the  material  contrasts  to  scattered  field  measurements 
taken  outside  the  object  domain.  Depending  on  the  inversion 
algorithm,  these  two  equations  are  used  in  combination  to 


recover  both  the  contrasts  and  the  total  fields.  Traditionally, 
( 1 )  and  (4)  are  used  as  they  are  to  develop  inverse  scattering 
algorithms. 

2.2.  Integral  Equations  for  Cavity  S- Parameter  Measurements. 
In  a  previous  work  [12],  we  showed  that  it  is  possible  to 
transform  (1)  and  (4)  so  that  they  are  consistent  with  an  S- 
parameter-based  measurement  system.  We  showed  that  the 
resulting  equations  were  valid  for  both  free-space  and  cavity¬ 
like  geometries  and  went  on  to  validate  the  free-space  case 
with  an  inverse  scattering  experiment  [13].  Here,  we  will 
summarize  the  results  for  a  cavity  geometry. 

Consider  the  cavity  depicted  in  Figure  1.  An  object  to 
be  imaged  is  placed  in  the  middle  of  the  cavity.  The  cavity 
is  filled  with  a  background  material  having  a  permittivity 
and  conductivity  of  £b  and  <7^,  respectively.  The  cavity  is 
lined  with  radiating  apertures,  which  could  be  antennas. 
Each  aperture  has  its  own  feeding  transmission  line  and  S- 
parameter  reference  plane. 
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We  define  the  normalized  incident  and  total  fields 
throughout  the  cavity  due  to  a  transmitting  aperture  as 


Cinc(1’) 

e(r) 


Ejnc(r) 

U0 

E(r) 

a0 


(6) 


where  a0  is  the  transmit  voltage  measured  with  respect  to 
the  5-parameter  reference  plane.  The  normalized  incident 
field  captures  all  background  multiple  scattering  not  present 
between  the  object  and  the  cavity. 

Let  transmitting  apertures  be  indexed  with  i  and  those 
receiving  indexed  with  j.  We  can  write  (1)  in  terms  of 
the  normalized  incident  and  total  fields  produced  by  a 
transmitter  by  dividing  both  sides  by  a0/. 


e,-( r)  =  einc,i(r)  +  J  G(r,r')  ■  0(r')e,(r')dV'.  (7) 

This  is  the  integral  equation  we  will  use  to  represent  the 
forward  scattering  solution.  The  normalized  total  field  is  the 
field  solution  in  the  object  domain  and,  with  the  appropriate 
dyadic  Greens  function  for  the  cavity,  includes  the  scattering 
interactions  between  the  object  and  the  cavity. 

In  [12]  we  showed  how  to  transform  the  scattered 
field  volume  integral  equation  given  by  (4)  into  one  that 
predicts  5-parameters.  This  new  integral  operator  allows 
us  to  directly  compare  model  predictions  to  measurements 
in  the  inversion  algorithm.  The  two-port  scattered  field  5- 
parameter,  SjiiSC a>  measured  between  the  transmission  line 
reference  planes  of  two  apertures  in  the  presence  of  an  object 
is  given  by 


SjUsa  =  J g j(r)  ■  0(r')ei(r')dV',  (8) 

where  ez(r)  is  the  normalized  total  object  field  produced  by 
the  transmitter  and  gj  (r)  is  the  vector  Green  s  function  kernel 
for  the  receiver.  It  was  also  shown  in  [12]  by  reciprocity  that 
g j  (r)  is  related  to  the  normalized  incident  field  of  the  receiver 
as 


g;«  =  - 


^inc,  / 
2lU)[l  J 


(r), 


(9) 


where  co  is  the  operating  frequency  in  radians,  [i  is  the  back¬ 
ground  permeability,  and  Z]0  is  the  characteristic  impedance 
of  the  receiver  transmission  line. 

Equations  (7)  and  (8)  are  the  integral  equations  we  will 
use  for  the  inverse  scattering  algorithm.  They  consistently 
link  the  electric  field  volume  integral  equations  to  an  5- 
parameter  measurement  system.  We  need  only  to  determine 
the  normalized  incident  fields  in  the  object  domain  and 
the  background  dyadic  Greens  function;  no  other  step  is 
required  to  characterize  the  system,  except  to  calibrate  the 
transmission  line  reference  planes. 

Lastly,  in  experiment,  we  never  measure  scattered  field 
5-parameters  directly  but  obtain  them  by  subtracting  the  5- 
parameters  for  the  total  and  incident  fields: 


S/z,sca  —  Sji,  tot  5jz>ino 


(10) 


Figure  9:  HFSS  model  of  a  simple  sphere  used  to  generate  synthetic 
scattered  field  S-parameters. 


where  5;  z,inc  is  measured  in  the  absence  of  the  object  and  5;  z>tot 
is  measured  in  the  presence  of  the  object. 


2.3.  Determining  einc(r)  and  G(r,r').  The  normalized  inci¬ 
dent  field  is  required  in  both  (7)  and  (9)  and  is  required 
for  every  aperture.  We  can  either  measure  it  experimentally 
or  estimate  it  with  simulation.  Experimentally  mapping  the 
fields  requires  proper  probe  calibration  and  has  the  added 
complication  in  a  cavity  that  the  probe-wall  interactions 
cannot  be  neglected.  An  alternative  approach,  the  one  we 
adopt  for  this  paper,  is  to  estimate  the  normalized  incident 
field  with  simulation.  This  can  be  done  provided  that  we 
have  a  computer  aided  design  (CAD)  model  that  accurately 
represents  the  cavity.  It  is  also  possible  in  simulation  to  model 
the  feeding  transmission  lines  and  line  voltages  in  order  to 
assign  an  5-parameter  reference  plane  that  is  identical  to 
the  reference  plane  used  by  a  vector  network  analyzer  for 
the  physical  measurement.  We  will  show  how  we  use  Ansoft 
FFFSS  to  accomplish  this. 

As  stated  in  the  introduction,  determining  the  back¬ 
ground  dyadic  Greens  function  is  nontrivial,  especially  for 
arbitrary  cavity  geometries.  Despite  this,  for  the  immediate 
investigation,  we  use  the  free-space  dyadic  Greens  func¬ 
tion  under  the  condition  that  the  background  medium  is 
extremely  lossy.  Though  not  strictly  correct,  this  approx¬ 
imation  is  convenient  provided  the  multiple  scattering 
throughout  the  cavity  is  limited  by  the  background  loss. 
It  also  allows  us,  for  the  time  being,  to  retain  use  of  an 
FFT-based  volumetric  forward  solver.  We  give  examples 
later  evaluating  this  assertion.  There  are  several  approaches 
for  determining  or  approximating  the  background  dyadic 
Greens  function  for  arbitrary  geometries,  which  we  discuss 
in  The  Appendix  and  leave  for  future  work. 
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BIM:  1 


Figure  10:  Reconstructions  of  a  single  sphere  (er,  o )  =  (40,  1)  located  at  (x,y,z)  =  (0,0,2  cm)  of  Example  1.  (a)  and  (c)  and  (b)  and  (d) 
Relative  permittivity  and  conductivity,  respectively,  (a)  and  (b)  is  the  Born  approximation,  (c)  and  (d)  is  BIM  iteration  4.  Here,  iterations 
help  retrieve  the  relative  permittivity  in  (c),  but  the  Born  approximation  yielded  better  conductivity  in  (b). 


3.  Born  Iterative  Method 

The  imaging  algorithm  we  use  is  the  Born  iterative  method 
(BIM)  [16-19] .  The  BIM  successively  linearizes  the  nonlinear 
problem  by  alternating  estimates  of  the  contrasts  and  the 
object  fields  according  to  the  following  algorithm. 

(1)  Assume  the  object  fields  are  the  incident  fields  (Born 
approximation). 

(2)  Given  the  measured  scattered  field  data,  estimate  the 
contrasts  with  the  current  object  fields  by  minimizing 
a  suitable  cost  function. 

(3)  Run  the  forward  solver  with  current  contrasts.  Store 
the  updated  object  field. 

(4)  Repeat  step  2  until  convergence. 

This  algorithm  and  its  implementation  are  described 
in  detail  in  our  previous  work  [13],  where  we  successfully 
formed  images  of  dielectric  constant  for  plastic  objects  in  a 
free-space  experiment.  This  was  done  using  the  same  BIM 


and  the  integral  equations  for  S-parameters  given  above  us¬ 
ing  antennas  characterized  with  HFSS. 

We  use  the  multivariate  covariance-based  cost  function 
of  [15].  The  Gaussian  interpretation  of  this  cost  function 
allows  us  to  experimentally  justify  the  values  we  use  to  regu¬ 
larize  it  by  our  a  priori  knowledge  of  the  experimental  noise 
and  range  of  contrast  values.  For  the  forward  solver,  because 
we  use  the  lossy  free-space  dyadic  Greens  function  to  model 
the  internal  scattering,  we  use  the  BCGFFT  [20-22],  which 
we  have  validated  with  analytic  solutions.  In  the  examples 
that  follow,  we  found  that  4  BIM  iterations  were  repeatably 
sufficient  for  the  data  residual  and  object  to  converge. 

4.  Breast  Imaging  System  Prototype 

The  breast  imaging  system  prototype  we  built  is  shown 
in  Figure  2.  The  imaging  structure  is  a  cavity,  shown  in 
Figure  3,  that  was  created  by  soldering  twelve  vertical  panels 
of  microwave  substrate  together  and  soldering  the  collection 
to  a  conducting  base.  Opposite  panels  are  separated  by  15  cm, 


International  Journal  of  Biomedical  Imaging 


7 


BIM:  1 


Figure  11:  Reconstructions  of  a  single  sphere  (er,  o )  =  (40,  0)  located  at  (x,y,z)  =  (0,0,2  cm)  of  Example  1.  Born  iterations  helped  retrieve 
the  relative  permittivity  in  (c)  and  are  essential  in  recovering  the  conductivity  in  (d). 


and  the  cavity  is  17  cm  long.  Three  antennas  are  printed  on 
each  panel  for  a  total  of  36  antennas.  In  the  prototype,  the 
three  antennas  of  one  panel  are  used  as  transmitters,  while 
all  other  antennas  are  receivers.  The  transmit  antennas  are 
switched  with  a  Dowkey  SP6T  electromechanical  switch.  The 
receivers  are  connected  through  an  SP33T  solid-state  switch¬ 
ing  matrix  that  was  designed  and  assembled  in-house.  2-port 
5-parameter  measurements  were  taken  with  an  Agilent  PNA- 
5230A  vector  network  analyzer  (VNA)  at  2.75  GHz  between 
each  transmitter  and  any  one  receiver.  This  frequency  was 
chosen  as  a  compromise  between  resolution  and  switch 
performance,  which  rolls  off  above  3  GHz.  A  rotator  was 
mounted  above  the  cavity  and  aligned  in  the  center  of  the 
cavity.  Test  objects  are  suspended  with  fishline  and  rotated  to 
provide  multiple  transmitter  views. 

4. 1 .  Liquid  Coupling  Medium.  We  expect  breast  tissue  to  have 
a  relative  permittivity  between  10  and  60  [23].  Without  a 
matching  medium,  much  of  the  incident  power  would  be 
reflected  at  the  breast/air  interface  reducing  the  sensitivity  of 
the  system  [24].  Also,  the  contrast  ratio  between  the  object 


and  the  background  would  be  too  high  for  the  BIM  inverse 
scattering  algorithm  to  converge. 

The  matching  medium  we  use  is  an  oil/ water  emulsion 
developed  in  a  previous  work  [25].  This  fluid  is  designed 
to  balance  the  high  permittivity  and  high  conductivity  of 
water  with  the  low  permittivity  and  low  conductivity  of  oil, 
in  order  to  achieve  a  fluid  with  moderate  permittivity  while 
limiting  loss  as  much  as  possible.  We  are  also  able  to  tune 
the  microwave  properties  of  this  emulsion  by  adjusting  the 
oil/water  ratio.  We  aimed  for  a  relative  permittivity  value 
around  20,  which  brings  the  maximum  permittivity  contrast 
to  about  3:1.  The  fluid  mixture  we  used  was  65%/35% 
oil/ water. 

The  electrical  properties  of  the  fluid  were  measured  us¬ 
ing  the  Agilent  85070E  slim  form  dielectric  probe.  The 
measured  properties  at  2.75  GHz  were  (er,<r)  =  (19,0.34). 
Relative  permittivity  is  unitless;  the  units  of  conductivity 
used  throughout  the  paper  are  Siemens/m.  When  using 
this  value  in  the  numerical  model  (presented  below)  the 
magnitude  of  cross -cavity  S21  required  some  adjustment 
when  compared  to  the  measurements.  We  obtained  the  best 
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BIM:  1 


Figure  12:  Reconstructions  of  a  single  sphere  (er,  a)  =  (10,  1)  located  at  (x,y,z)  =  (0,0,2  cm)  of  Example  1.  Born  iterations  helped  the 
recovery  of  the  low  permittivity  in  (c),  but  at  the  expense  of  the  correct  conductivity  value  which  was  better  with  the  Born  approximation  in 

(b). 


model  agreement  for  (er,a)  =  (21,0.475),  which  are  the 
values  we  use  throughout  the  paper.  We  suspect  that  the 
probe  area  may  be  too  small  to  accurately  measure  the  bulk 
properties  of  the  mixture,  but  the  fluid  otherwise  appears 
homogeneous  for  propagation  at  2.75  GHz.  We  are  still 
investigating  this  effect. 

When  taking  data,  we  fill  the  cavity  with  the  coupling 
fluid  to  a  height  that  is  0.5  cm  below  the  top  edge.  This 
fluid  height  is  accounted  for  in  the  numerical  model.  Any 
fluid  displacement  from  adding  or  removing  test  objects 
is  compensated  in  order  to  keep  the  height  constant. 
We  have  also  found  the  emulsion  to  be  stable  over  the 
course  of  measurements,  which  we  confirmed  by  comparing 
transmission  measurements  before  and  after  we  take  data  for 
imaging. 

4.2.  Antenna  Design.  The  antennas  are  bow-tie  patch  anten¬ 
nas,  similar  to  the  antennas  in  [6,  26].  They  are  of  single 
frequency  and  vertical  polarization.  The  bow-tie  antenna 
was  chosen  to  give  more  degrees  of  freedom  to  help 


impedance  match  the  antenna  to  the  coupling  fluid.  The 
vertical  polarization  was  chosen  for  best  illumination  of  the 
object  and  other  antennas  in  the  cylindrical  geometry.  The 
substrate  material  is  Rogers  RO3210,  with  50  mil  thickness 
and  reported  dielectric  constant  of  10.2.  The  antennas  were 
originally  designed  to  operate  at  2.8  GHz  in  the  cavity  filled 
with  a  fluid  with  (er,  a)  =  (24, 0.34);  however,  after  iterating, 
we  found  best  performance  at  2.75  GHz  in  a  fluid  of  (er,  o )  = 
(21,0.475). 

4.3.  System  Parameters.  In  determining  the  system  noise 
and  isolation  requirements,  the  minimum  expected  signal 
determines  the  required  noise  level,  and  the  maximum 
relative  magnitude  between  signals  on  adjacent  channels 
determines  the  required  switch  path  isolation.  From  previous 
numerical  studies  of  cavity-like  breast  imaging  with  similar 
emulsion  properties  [27],  we  expect  the  scattered  field  S21 
magnitude  of  small  inclusions  to  be  in  the  range  from 
-100  to  -50  dB,  and  so  the  relative  signal  strength  between 
adjacent  antennas  could  differ  by  as  much  as  -50  dB.  This 
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BIM:  1 


Figure  13:  Reconstructions  of  a  single  sphere  (er,  o)  =  (40,  0)  located  at  (x,y,z)  =  (0,0,2  cm)  of  Example  1.  Born  iterations  helped  bring 
out  the  proper  conductivity  value  in  (d). 


means  that  the  noise  of  our  system  must  be  less  than 
-lOOdB,  which  is  achievable  by  our  VNA  with  averaging 
and  an  IF  bandwidth  of  100  kHz  or  less.  Also,  the  switching 
matrix  paths  must  be  isolated  by  at  least  -50  dB. 

4.4.  Switching  Matrix.  The  receivers  were  connected  through 
a  SP33T  solid-state  switching  matrix  that  was  designed  and 
assembled  in-house.  The  matrix  consists  of  two  custom 
SP16T  solid-state  switching  matrices  and  a  cascaded  pair  of 
Miniciruits  SPDT  switches.  Each  SP16T  switch  is  composed 
of  two  layers  of  SP4T  Hittite  HMC241QS16  nonreflective 
switches,  which  are  buffered  at  the  output  by  a  third  layer 
consisting  of  a  single  SPDT  Hittite  HMC284MS8GE  on 
each  path.  The  buffer  layer  was  added  to  increase  interpath 
isolation.  The  switch  is  controlled  with  an  embedded  digital 
board  and  computer  parallel  port.  The  operating  band  of 
the  switching  matrix  is  between  0.1-3  GHz.  The  overall  loss 
of  a  path  through  the  SP33T  matrix  is  no  worse  than  8  dB 
across  the  band.  We  measured  the  switch  path  isolation  to 
be  better  than  -55  dB  between  1-3  GHz,  which  meets  the 
criteria  above. 


By  separating  the  transmitter  and  receiver  switching,  the 
isolation  between  these  two  operation  modes  is  dictated 
by  the  network  analyzer  and  the  cables.  In  more  realistic 
systems,  where  the  antennas  are  dual  mode  and  so  object 
rotation  is  not  necessary,  the  isolation  requirements  are  more 
stringent,  because  the  transmit  amplitude  will  be  orders  of 
magnitude  larger  than  the  scattered  field. 

4.5.  VNA  Calibration.  Two-port  VNA  calibrations  were 
accomplished  between  each  transmitter  and  each  receiver. 
The  S-parameter  reference  planes  were  calibrated  to  the 
points  where  the  cables  connect  to  the  antenna.  These 
reference  planes  are  identical  to  those  in  the  HFSS  CAD 
model  (presented  below).  While  calibrating,  we  left  the 
unused  ports  open  with  the  rationale  that  the  one-way 
switch  isolation  of  -55  dB  provided  sufficient  matching  to 
the  open  ports.  Short-open-load  measurements  for  a  1- 
port  calibration  were  taken  for  each  antenna.  Next,  we 
measured  the  through  path  between  the  transmitter  and  each 
receiver  using  a  connector.  In  software,  we  combined  the 
1-port  and  through  measurements  to  accomplish  a  2-port 
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short- open-load- through  (SOLT)  calibration  with  arbitrary 
through  between  each  transmitter  and  receiver  for  a  total 
of  99  separate  2 -port  calibrations.  The  calibration  for  a 
particular  transmitter/receiver  pair  is  recalled  in  the  VNA 
before  taking  data. 

5.  Numerical  Model 

We  use  Ansoft  HFSS  to  numerically  model  the  cavity,  similar 
to  [27].  We  use  it  several  ways.  First,  we  model  the  feeding 
transmission  lines  in  order  to  assign  5-parameter  reference 
planes  that  are  identical  in  both  simulation  and  experiment. 
Second,  we  estimate  the  normalized  incident  fields  due  to 
the  transmitters  throughout  the  cavity  for  use  in  (7)  and 
(9),  where  the  normalized  incident  fields  now  include  all 
background  multiple  scattering  not  present  between  the 
object  and  the  cavity.  Also,  we  use  the  model  to  generate 
synthetic  scattered  field  S-parameters  of  numerical  targets 
in  order  to  study  the  performance  of  the  inverse  scattering 
algorithm  given  the  source  geometry  and  system  parameters. 

Figure  4  shows  the  HFSS  CAD  model  of  the  12-sided 
cavity.  The  model  includes  the  panel  thickness  and  dielectric 
constant,  bottom  conductor,  probe  feed,  coupling  fluid  prop¬ 
erties,  and  height  of  the  fluid.  Same  as  in  the  experiment,  the 
cavity  is  filled  to  a  height  that  is  0.5  cm  below  the  top  (seen 
as  the  line  below  the  top  edge  of  the  cavity).  The  remaining 
0.5  cm  is  air  with  a  radiating  boundary  condition.  The  outer 
boundary  of  the  cavity  is  PEC. 

Next  we  compare  measured  and  simulated  incident  S- 
parameters  in  order  to  access  the  accuracy  of  the  model. 
Figure  5  shows  the  magnitude  and  phase,  respectively,  of 
the  measured  and  simulated  incident  S-parameters  between 
each  transmitter  and  all  receivers.  The  magnitude  and  phase 
agree  best  when  the  receivers  are  on  the  same  level  as  the 
transmitter.  In  this  case,  the  magnitude  agrees  generally  to 
within  3  dB,  for  all  three  levels,  and  the  phase  agrees  to  within 
30  degrees,  which  is  approximately  A/ 10,  a  common  metric 
for  many  microwave  systems.  For  measurements  between 
antenna  levels  in  Figure  5,  the  agreement  is  not  as  good 
in  magnitude,  but  the  phase  error  remains  similar  to  the 
previous  cases.  This  also  shows  that  the  one-way  path  loss 
across  the  cavity  is  approximately  -50  dB,  so  we  expect  any 
multiple  scattering  to  be  localized.  This  partially  justifies  our 
approximation  of  the  cavity  dyadic  Green  s  function  with  the 
lossy  free-space  dyadic  Greens  function. 

When  computing  the  incident  fields,  the  center  of  the 
cavity  was  meshed  with  a  coarse  Cartesian  grid  of  sparse 
unassigned  sheets,  shown  in  Figure  6.  Sheets  are  spaced 
every  5  mm  in  the  x,  y,  and  z  directions.  The  spacing  is 
approximately  A/5  at  2.75  GHz  in  the  fluid  with  relative 
permittivity  of  21.  We  have  found  that  this  helps  constrain 
the  adaptive  meshing  of  HFSS  when  we  obtain  the  incident 
fields  by  interpolating  the  FEM  mesh  onto  a  fine  Cartesian 
grid,  [12]. 

When  simulating  the  structure,  with  or  without  scat¬ 
tering  targets,  we  use  a  convergence  criterion  of  AS  = 
0.02  which  is  reached  in  7  adaptive  meshing  iterations. 
A  typical  simulation  completed  with  approximately  1.4 
million  tetrahedra  using  23.5  GBytes  of  RAM  and  swap 
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Figure  14:  HFSS  CAD  numerical  breast  phantom  of  Example  2. 
The  inclusion  is  2  cm  in  diameter  with  relative  permittivity  and 
conductivity  contrasts  of  2  : 1 .  The  skin  layer  is  2  mm  thick. 


space  to  obtain  a  full  36  x  36  S-matrix.  Simulations  took 
approximately  25  hours  on  a  dual  E5504  Intel  Xeon  (2x  Quad 
Core)  desktop  with  24  GBytes  of  RAM.  Figure  7  shows  a 
typical  convergence  rate  as  a  function  of  tetrahedra. 

We  obtained  the  incident  fields  for  only  the  three 
transmitters.  The  incident  fields  for  the  receivers  were 
obtained  through  rotation,  where  we  assume  the  12  panels 
of  the  experimental  cavity  are  identical.  The  incident  fields 
were  sampled  on  a  17  cm  X  17  cm  x  18  cm  grid  with  1  mm 
spacing,  which  is  A/24  at  2.75  GHz  in  a  fluid  with  relative 
permittivity  21.  In  simulation,  the  average  transmit  power 
was  1  Watt,  so,  from  transmission  line  analysis,  the  line 
voltage  is  given  by 

a0  =  sjlP^Zo  =  'fiz'o,  (11) 

which  is  used  in  (6).  The  phase  of  a0  is  zero  because  the 
5-parameter  reference  planes  of  the  HFSS  model  and  the 
experimental  cavity  were  identical. 

Figure  8  shows  three  crosscuts  of  the  z- component  of  the 
incident  electric  field  through  the  center  of  the  cavity  for  the 
center  transmitter  in  a  fluid  of  relative  permittivity  of  21  and 
conductivity  0.475  at  2.75  GHz.  The  coordinate  origin  is  at 
the  center  of  the  cavity,  and  the  transmitter  is  located  on  the 
positive  x  axis.  The  effects  of  the  cavity  on  the  incident  field 
are  seen  in  Figure  7,  where  the  fields  are  guided  by  the  walls 
of  the  cavity;  the  coaxial  feeds  are  also  visible;  the  fluid- air 
interface  is  visible  in  Figures  8(b)  and  8(c). 

6.  Image  Reconstructions 

6.1.  Synthetic  Data.  We  first  test  the  BIM  and  numerical 
characterization  using  synthetic  data  from  HFSS.  This  is 
to  assess  the  performance  of  the  algorithm  and  source 
geometry  under  near  ideal  circumstances.  We  simulated  the 
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Figure  15:  Reconstructions  of  the  HFSS  numerical  breast  phantom  in  Example  2  after  four  iterations,  (a),  (c),  and  (e)  and  (b),  (d),  and  (f) 
Relative  permittivity  and  conductivity,  respectively,  (a)  and  (b),  (c)  and  (d),  and  (e)  and  (f)  Cuts  at  x  =  0  cm,  y  -  0  cm,  and  z  =  3  cm.  The 
relative  permittivity  of  the  inclusion  is  recovered,  but  both  images  contain  many  artifacts. 
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Figure  16:  HFSS  CAD  numerical  breast  phantom  with  skin 
layer,  fat  layer,  glandular  tissue,  and  chest  wall  of  Example  3. 
The  inclusion  is  1  cm  in  diameter  with  relative  permittivity  and 
conductivity  contrasts  of  2  : 1 . 


scattered  field  S-parameters  of  simple  numerical  objects  and 
use  these  data  as  measurements  in  the  inversion  algorithm. 
FFFSS  scattered  field  data  includes  any  multiple  scattering 
between  the  object  and  the  cavity.  The  background  medium 
had  a  relative  permittivity  of  21  and  a  conductivity  of 
0.475  Siemens/m.  The  incident  fields  were  computed  with 
these  background  parameters  and  used  in  volume  integral 
equations. 

Example  1.  We  first  used  FFFSS  to  simulate  the  scattered 
field  S-parameters  for  a  single  1.5  cm  diameter  sphere  located 
at  (x,y,z)  =  (0,0,2  cm)  with  four  combinations  of  relative 
permittivity  and  conductivity:  (40,  1),  (40,  0),  (10,  1),  and 
(10,  0).  The  FFFSS  model  is  shown  in  Figure  9.  Figures  10, 
11,  12,  and  13  show  images  of  the  first  and  fourth  BIM 
iterations  for  each  object.  As  shown,  in  some  cases,  the 
BIM  steps  were  essential  in  recovering  the  correct  property 
values  of  the  sphere;  in  other  cases,  the  relative  permittivity 
was  improved  at  the  expense  of  the  conductivity  value. 
These  images  show  that  the  source  geometry  and  numerical 
characterization  are  adequate  for  the  retrieval  of  some  object 
property  combinations,  but  not  others.  This  fact,  together 
with  the  visible  artifacts,  suggests  that  the  images  could  be 
improved  with  a  denser  source  geometry. 

Example  2.  Next  we  imaged  a  more  anatomical  numerical 
breast  phantom.  The  numerical  phantom  is  shown  in 
Figure  14.  The  breast  is  9  cm  at  the  widest  point  and  6  cm 
deep.  The  outer  layer  is  a  2  mm  thick  skin  layer,  and  the 
inclusion  is  2  cm  in  diameter.  The  dielectric  properties  of  the 
skin  layer,  glandular  tissue,  and  inclusion,  respectively,  are 
(er,  o)  =  {(45,1.59),  (21,  0.475),  and  (42,  0.8)},  which  were 
obtained  from  [28].  We  assume  we  know  the  volume  region 
of  the  breast,  so  we  mask  that  volume  excluding  all  other 


(b) 


Figure  17:  Reconstructions  of  the  HFSS  numerical  breast  phantom, 
which  includes  the  chest  wall  of  Example  3.  (a)  and  (b)  Relative 
permittivity  and  conductivity,  respectively.  The  object  could  not  be 
reconstructed. 


points  during  inversion.  Figure  15  shows  the  reconstructed 
relative  permittivity  and  conductivity  after  4  iterations  for 
three  cuts.  The  relative  permittivity  of  the  inclusion  is 
recovered,  but  the  conductivity  of  the  inclusion  is  not 
recovered.  The  skin  layer  is  also  visible  in  the  conductivity 
images.  Both  sets  of  images  suffer  from  artifacts,  which  is  due 
to  the  sparse  spatial  sampling  of  the  antennas  and  indicates 
that  the  images  can  be  improved  with  more  angular  views. 

Example  3.  To  push  the  algorithm,  we  imaged  a  phantom 
that  included  a  skin  layer,  fat  layer,  glandular  tissue,  chest 
wall,  and  inclusion,  with  relative  permittivity  and  conductiv¬ 
ity,  respectively,  of  (45,  1.6),  (5.1,0.16),  (21,0.475),  (52,  2.0), 
and  (40,  1.0).  The  HFSS  model  is  shown  in  Figure  16.  The 
reconstructions  are  shown  in  Figure  17.  In  this  case  the 
algorithm  failed  to  recover  the  contrasts.  This  suggests  that 
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Figure  18:  Continued. 
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BIM:  4  eb  =  22 


(b) 


Figure  18:  Sensitivity  of  image  reconstructions  to  background  permittivity  for  Example  4.  (a),  (c),  (e),  (g),  and  (i)  and  (b),  (d),  (f),  (h),  and 
(j)  Relative  permittivity  and  conductivity,  respectively.  The  scattered  field  data  was  generated  in  HFSS  in  a  background  of  (21,  0.475).  The 
reconstructions  are  done  with  assumed  background  relative  permittivities  of  20,  20.5,  21,  21.5,  22  for  (a)  and  (b),  (c)  and  (d),  (e)  and  (f), 
(g)  and  (h)  and  (i)  and  (j),  respectively.  The  recovered  contrasts  of  the  sphere  oscillate  about  the  background. 


(1)  the  object  is  too  different  from  the  background  for  the 
BIM  to  converge,  (2)  object-cavity  interactions  are  too  strong 
to  use  the  free- space  dyadic  Greens  function,  or  (3)  images 
cannot  be  constructed  if  the  chest  wall  is  not  modeled, 
meaning  that  it  is  necessary  to  model  the  chest  wall  for  the 
incident  fields  and  the  dyadic  Green  s  function. 

Example  4.  Finally,  we  studied  the  effects  of  different  back¬ 
ground  permittivities  when  forming  images.  This  represents 
a  case  in  experiment  where  the  measurements  are  taken  in 
a  fluid  with  some  set  of  properties,  but  the  fluid  properties 
we  use  in  the  model  are  slightly  off.  We  formed  images  using 
FFFSS  scattered  field  data  of  the  sphere  with  (er,  o)  =  (40,  0) 
in  a  background  of  (e&,  o)  =  (21,  0.475),  but  where  we  use 
incident  fields  from  five  different  background  permittivities: 
{20,  20.5,  21  (again),  21.5,  22}  and  the  same  conductivity. 

Figure  18  shows  3D  crosscuts  at  the  fourth  BIM  iteration 
for  all  five  backgrounds.  Figures  18(e)  and  18(f)  are  the 
correct  images.  Notice  that  an  error  in  the  background 
permittivity  of  1,  or  5%,  is  enough  for  the  reconstructed 
object  contrast  to  oscillate,  demonstrating  that  reconstruc¬ 
tions  are  very  sensitive  to  our  knowledge  of  the  background 
properties. 

6.2.  Experimental  Data.  At  this  time,  only  simple  plastic 
objects  have  been  imaged  with  the  experimental  system; 
however,  future  work  includes  imaging  more  realistic  breast 
phantoms.  Among  the  test  objects,  we  show  the  results 
here  for  several  acrylic  spheres.  The  objects  were  suspended 
from  a  platform  and  rotated  to  12  positions  in  30  degree 
increments.  Scattered  field  S-parameter  measurements  from 
each  position  were  combined  to  yield  a  full  36  X  36  S- 
parameter  matrix,  which  was  used  in  the  inverse  scattering 
algorithm. 

Experiment  1.  We  imaged  a  single  acrylic  sphere,  shown  in 
Figure  19.  The  diameter  of  the  sphere  was  2.54  cm,  with 


properties  (er,  a)  =  (2.7,  0).  The  sphere  was  located  at 
approximately  (x,  y ,  z)  =  (1.5  cm,  1.5  cm,  0).  Figure  20  shows 
the  reconstructions  after  4  iterations  of  the  x-y  plane. 
The  inversion  domain  is  masked  so  that  only  a  cylindrical 
region  containing  the  rotated  object  is  imaged.  We  also 
imaged  two  acrylic  spheres,  shown  in  Figure  19.  Figure  21 
shows  the  reconstructions  after  4  iterations.  In  both  cases, 
the  relative  permittivity  is  recovered  quite  well,  and  the 
conductivity  contrast  is  correctly  valued  but  the  shape  is 
incorrect.  There  are  also  many  artifacts  present.  Given  that 
the  imaging  algorithm  could  recover  the  single  sphere  using 
FFFSS  data,  we  can  attribute  these  discrepancies  to  differences 
between  the  experiment  and  the  model,  such  as  knowledge  in 
the  coupling  medium  properties,  substrate  properties,  VNA 
calibration,  cavity  size  measurements,  or  object  motion. 

Experiment  2.  Finally,  while  the  primary  discussions  in  this 
paper  concern  a  cavity  having  antennas  that  operate  at 
2.75  GFFz,  we  also  built  a  lower  frequency  cavity  where  the 
antennas  operate  at  915MFFz.  This  cavity  was  numerically 
characterized  using  the  same  methods,  but  the  background 
fluid  properties  were  (er,  a)  =  (23,  0.1).  Figure  22  shows 
the  cavity  with  three  acrylic  spheres.  Two  spheres  are 
located  in  the  x-y  plane,  while  the  third  is  positioned  at 
approximately  (x,y,z)  =  (4  cm,  -3  cm,  5  cm).  We  imaged 
the  relative  permittivity  and  conductivity,  and  the  results 
after  4  iterations  are  shown  in  Figure  23.  The  shape  and 
properties  of  the  two  in-plane  spheres  are  well  recovered. 
The  third  sphere  is  also  detected  but  cut  off  at  the  upper 
left  of  the  imaging  domain.  Artifacts  are  also  present,  but 
this  example  better  demonstrates  that  the  numerical  charac¬ 
terization,  BIM,  and  free-space  Green  s  function  are  capable 
of  recovering  objects  in  this  cavity  and  source  geometry.  It 
should  be  noted  that  images  formed  with  data  at  915MFFz 
are  less  susceptible  to  modeling  errors  because  the  cavity 
and  objects  are  electrically  smaller,  but  the  resolution  is 
reduced. 
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Figure  19:  Test  objects  and  coupling  fluid  for  Experiment  1.  (a) 
Single  suspended  acrylic  sphere,  (b)  Two  acrylic  spheres,  (c)  Cavity 
filled  with  the  coupling  medium.  Objects  are  suspended  and  rotated 
from  the  nylon  platform. 


6.3.  Discussion.  Overall,  the  imaging  algorithm,  numerical 
characterization,  and  experiment  worked  with  some  success, 
and  there  are  several  areas  for  continued  investigation. 

First,  Examples  1  and  2,  and  also  Experiments  1  and 
2,  validate  the  technique  described  in  this  paper  showing 
that  the  numerical  characterization  of  the  cavity  incident 
fields  and  the  use  of  the  vector  Green  s  function  formulation 
linking  the  incident  fields  to  the  inverse  scattering  algorithm 
can  be  used  to  successfully  form  images  in  a  cavity  geometry. 
Examples  1  and  2  demonstrate  the  consistency  of  the  method 
using  synthetic  scattered  field  5-parameter  data.  Experiments 
1  and  2  show  that  the  characterization  and  experiment 
agreed  enough  for  the  BIM  to  recover  the  location  and 
permittivity  of  the  test  objects.  More  realistic  phantoms 
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Figure  20:  Reconstructions  of  the  single  acrylic  sphere  of 
Experiment  1  shown  in  Figure  19(a).  (a):  Relative  permittivity,  (b): 
Conductivity.  The  permittivity  is  recovered  well  but  the  shape  in  the 
conductivity  is  not. 


and  lower  contrast  phantoms  will  help  further  confirm  the 
methodology. 

Second,  in  Example  1,  although  some  permittivity  and 
conductivity  combinations  of  the  sphere  were  recovered, 
others  were  not.  Given  that  the  data  was  synthetic,  this 
points  to  inherent  imaging  ambiguities  in  the  simultaneous 
retrieval  of  both  permittivity  and  conductivity  in  the  inverse 
scattering  problem.  Possible  solutions  are  increasing  the 
number  of  unique  data,  or  including  prior  information 
about  the  relations  between  permittivity  and  conductivity  in 
tissue. 

Third,  the  success  of  the  algorithm  in  Example  2  in 
recovering  the  partial  breast  phantom  suggests  that  our  use 
of  the  lossy  free-space  dyadic  Greens  function  in  the  forward 
solver  of  the  BIM  did  not  grossly  affect  image  reconstruction 
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Figure  21:  Reconstructions  of  the  two  acrylic  spheres  of 
Experiment  1  shown  in  Figure  19(b).  (a)  Relative  permittivity,  (b) 
Conductivity.  The  permittivity  is  again  recovered  well  but  the  shape 
in  the  conductivity  is  not. 


in  this  case.  This  is  keeping  in  mind  that  the  synthetic 
scattered  field  5-parameter  data  did  include  any  multiple 
scattering  between  the  phantom  and  the  cavity. 

Fourth,  in  light  of  the  successful  reconstruction  of  the 
simple  phantom  in  Example  2,  the  failure  of  the  algorithm 
to  recover  the  more  complete  breast  phantom  in  Example  3 
points  to  the  need  to  model  the  chest  wall.  This  can  be 
done  by  including  it  in  the  incident  field  computations  but  it 
may  also  be  necessary  in  estimating  the  cavity  dyadic  Greens 
function.  This  is  an  area  to  be  investigated. 

Lastly,  Example  4  shows  that  we  must  know  the  back¬ 
ground  relative  permittivity  to  within  5%  of  the  actual  or  else 
risk  incorrectly  estimating  whether  the  contrasts  are  higher 
or  lower  than  the  background.  An  equivalent  error  can 
arise  from  a  correct  background  permittivity  but  incorrectly 


Figure  22:  Second  cavity  with  antennas  designed  to  operate  at 
915  MHz  of  Experiment  2.  Three  acrylic  sphere  are  suspended  (one 
visible).  Cavity  is  filled  with  fluid  for  imaging. 


measuring  the  dimensions  of  the  cavity.  We  suspect  that  the 
very  high  recovered  conductivity  values  in  both  Experiments 
1  and  2  may  be  due  in  part  to  these  types  of  systematic  errors. 
This  demonstrates  the  difficulty  in  achieving  the  necessary 
consistency  between  the  model,  experiment,  characteriza¬ 
tion,  and  imaging  algorithm  to  accurately  form  microwave 
breast  images  of  diagnostic  quality. 

7.  Conclusion 

We  demonstrated  the  use  of  a  numerical  characterization 
technique  for  a  breast  imaging  system  prototype.  We 
used  HFSS  to  numerically  estimate  the  incident  fields  of 
the  antennas  in  a  cavity  geometry  and  formally  linked 
them  to  an  S-parameter-based  inverse  scattering  algorithm 
and  experimental  setup.  The  imaging  algorithm  was  the 
Born  Iterative  Method  and  recovered  both  numerical  and 
experimental  test  objects  with  some  success.  Future  work 
includes  further  validation  of  our  methodology,  imaging 
realistic  breast  phantoms,  investigating  practical  solutions  to 
modeling  breast- cavity  scattering  interactions,  image  quality 
assessments  with  and  without  numerical  characterization, 
and  developing  a  hemispherical  cavity  and  clinical  imaging 
system. 

Appendix 

Determining  the  Background  Dyadic 
Green's  Function 

We  list  the  following  approaches  for  obtaining  the  back¬ 
ground  dyadic  Greens  function  as  work  for  future  investi¬ 
gation. 

Analytical  Dyadic  Greens  Function.  There  exist  analytical 
solutions  of  the  dyadic  Green  s  function  for  some  simple  cav¬ 
ity  geometries,  such  as  cubes  or  cylinders,  [29],  which  might 
approximately  model  certain  cavity-based  imaging  setups. 
These  solutions,  however,  will  likely  not  include  finer  details 
such  as  antenna  plating,  connectors,  substrate  material,  or 
open-ended  cavities,  such  as  those  used  for  breast  imaging. 
Analytic  solutions  though  lend  themselves  to  the  possibility 
of  retaining  some  convolution  structure  in  the  VIE  so  fast 
forward  solvers  can  be  used  (e.g.,  fast  half  space  solutions 
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Figure  23:  Reconstruction  of  the  relative  permittivity  and  conductivity  of  three  acrylic  spheres  using  a  cavity  operating  at  915  MHz  from 
Experiment  2.  The  two  spheres  in  plane  are  well  recovered  and  the  third  detected  at  the  upper  left  of  the  image. 


[30],  applied  to  multisided  cavities).  Determining  the  dyadic 
Greens  function  analytically  becomes  a  formidable  task  as 
the  geometry  complicates,  where  simulation  may  be  better 
suited. 

Full  Numeric  Simulation.  The  most  complete  solution 
is  to  fully  simulate  the  object  and  cavity  using  a  numeric 
simulator,  which  will  capture  all  the  multiple  scattering 
between  the  object  and  the  cavity.  However,  unlike  a  dyadic 
Green  s  function,  which  only  needs  to  be  found  once  for  a 
particular  geometry  and  the  values  of  which  are  only  needed 
on  the  interior  of  the  object  domain,  this  method  must 
simulate  the  cavity  structure  outside  the  object  domain  in 
every  instance  of  the  simulation.  When  used  in  an  inverse 
scattering  algorithm,  which  might  compute  the  domain  VIE 
for  each  source,  frequency,  and  iteration,  then  repeatedly 
simulating  the  cavity  structure  adds  to  the  already  high 
computational  burden.  In  addition,  one  must  choose  a 
proper  simulation  technique  to  handle  both  antenna  surfaces 
and  inhomogenous  media. 

Numerical  Dyadic  Greens  Function.  If  analytical  solutions 
are  not  accurate  enough,  then  one  must  determine  the  dyadic 
Green  s  function  numerically.  This  requires  simulating  three 
orthogonal  dipoles  in  turn  at  every  point  in  the  object 
domain  and  recording  the  response  at  every  other  point 
in  the  domain.  The  dyadic  Greens  function  is  symmetric, 
so  half  of  the  combinations  are  redundant,  and  while  the 
convolution  nature  of  the  VIE  is  destroyed,  some  compu¬ 
tational  speed-up  is  possible  for  symmetric  operators.  This 
technique,  however,  requires  accurate  modeling  around  the 
dipole  singularity,  which  can  be  difficult.  In  the  case  of  PEC 
structures,  the  technique  in  [31]  computes  the  dyadic  Greens 
function  by  finding  an  array  of  image  dipoles  outside  the 
cavity,  which  avoids  the  complications  from  the  singularity. 
The  main  advantage  of  determining  the  dyadic  Greens 
function  numerically  is  that,  once  found,  we  no  longer  need 
to  simulate  the  cavity  structure  and  can  turn  our  attention  to 
optimizing  the  computation  of  the  dyadic  Greens  function. 


Approximate  Solutions.  If  the  background  loss  is  suffi¬ 
ciently  high,  so  that  the  resonances  of  the  cavity  are  damped, 
then  we  can  approximate  the  dyadic  Green’s  function.  This 
can  be  done  by  adopting  an  analytical  solution  (e.g.,  free- 
space  or  cavity)  or  by,  for  instance,  developing  a  perturbation 
method.  Adopting  the  free-space  dyadic  Green’s  function  (or 
a  perturbation  on  it)  also  allows  us  to  retain  the  convolution 
structure  of  the  VIE-  and  FFT-based  forward  solvers,  which 
may  be  more  beneficial  to  the  inverse  scattering  algorithm 
than  modeling  higher-order  multiple  scattering. 
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Abstract — A  self-contained  source  characterization  method 
for  commercial  ultrasound  probes  in  transmission  acoustic 
inverse  scattering  is  derived  and  experimentally  tested.  The 
method  is  based  on  modified  scattered  field  volume  integral 
equations  that  are  linked  to  the  source-scattering  transducer 
model.  The  source-scattering  parameters  are  estimated  via 
pair-wise  transducer  measurements  and  the  nonlinear  inver¬ 
sion  of  an  acoustic  propagation  model  that  is  derived.  This 
combination  creates  a  formal  link  between  the  transducer  char¬ 
acterization  and  the  inverse  scattering  algorithm.  The  method 
is  tested  with  two  commercial  ultrasound  probes  in  a  transmis¬ 
sion  geometry  including  provisions  for  estimating  the  probe 
locations  and  aligning  a  robotic  rotator.  The  transducer  char¬ 
acterization  results  show  that  the  nonlinear  inversion  fit  the 
measured  data  well.  The  transducer  calibration  and  inverse 
scattering  algorithm  are  tested  on  simple  targets.  Initial  imag¬ 
es  show  that  the  recovered  contrasts  are  physically  consistent 
with  expected  values. 


I.  Introduction 

Inverse  scattering  has  received  attention  as  a  means  of 
improving  medical  ultrasound  by  providing  quantita¬ 
tive  images  of  tissue  material  properties,  such  as  speed 
of  sound,  without  the  traditional  artifacts  of  backscatter 
ultrasound  [1],  [2].  Inverse  scattering  systems,  however, 
require  full-wave  scattering  simulations,  large-scale  op¬ 
timization  algorithms,  and  absolute  source  characteriza¬ 
tion.  Despite  ongoing  progress  in  system  and  algorithm 
development  [3]— [10] ,  transducer  characterization  remains 
an  open  topic  in  making  ultrasound  inverse  scattering  ex¬ 
perimentally  viable.  The  effects  of  incomplete  transducer 
characterization  on  image  quality  in  full-wave  breast  ul¬ 
trasound,  for  example,  have  been  reported  in  [11]. 
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Not  all  transducer  characterization  methods  are  suit¬ 
able  for  full-wave  imaging.  Phased-array  beamforming,  for 
instance,  requires  only  magnitude  radiation  patterns  mea¬ 
sured  with  hydrophones.  Plane  wave  approximations  are 
also  sufficient  for  many  geometries.  Also,  recent  demon¬ 
strations  of  hybrid  speed  of  sound  and  diffraction  tomog¬ 
raphy  algorithms  have  successfully  produced  quantitative 
ultrasound  images  without  transducer  characterization 
[10],  [12].  Full- wave  inverse  scattering,  however,  requires 
transducer  field  characterization  with  absolute  phase;  the 
transducer  model  must  relate  the  actual  radiated  fields  to 
the  input / output  voltages  using  more  than  an  ideal  trans¬ 
fer  function.  Even  if  a  full-wave  characterization  exists, 
it  is  not  always  clear  how  to  incorporate  the  information 
into  the  imaging  algorithms.  Thus,  we  seek  a  transducer 
characterization  method  that  can  capture  the  full-wave 
behavior  of  a  transducer  and  can  also  be  rigorously  linked 
to  the  equations  of  the  inverse  scattering  algorithm. 

In  previous  work,  we  developed  a  formalism  for  micro- 
wave  inverse  scattering  and  antenna  characterization  in 
which  we  linked  the  source  model,  inversion  algorithm,  ob¬ 
ject  properties,  and  measurements  in  a  consistent  frame¬ 
work  [13],  [14].  We  proved  the  method  in  several  experi¬ 
ments  [15],  [16]  in  which  we  formed  quantitative  images  of 
dielectric  constant  without  the  use  of  calibration  targets. 
Given  the  success  of  this  method  in  microwave  imaging, 
we  sought  to  develop  and  test  the  analogous  method  for 
ultrasound  inverse  scattering. 

Although  the  derivation  of  the  acoustic  formalism  is 
straightforward,  key  differences  between  microwave  and 
ultrasound  experimentation  have  motivated  several  ad¬ 
ditional  developments.  First,  much  of  the  success  in  our 
microwave  experiments  was  due  to  our  ability  to  fully  sim¬ 
ulate  antenna  sources  using  commercial  electromagnetic 
software  packages.  This  allowed  antenna  characterization 
to  be  independent  of  the  measurement  setup.  We  are  un¬ 
able  to  do  this  with  ultrasonic  transducers  given  their 
small  size  and  unknown  construction.  Without  full  elec¬ 
tromechanical  simulation,  we  must  characterize  the  trans¬ 
ducers  from  measurement  alone.  Second,  performing  hy¬ 
drophone  measurements  that  satisfy  the  requirements  for 
full-wave  characterization  is  not  trivial  and  is  complicated 
by  the  effects  of  the  hydrophone.  Assuming  that  simula¬ 
tion  and  hydrophone  measurements  are  unavailable,  we 
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would  still  like  to  characterize  the  transducers  within  the 
framework  we  are  going  to  develop. 

These  considerations  and  our  current  interest  in  trans¬ 
mission  imaging  have  prompted  us  to  explore  a  new  self¬ 
characterization  technique  for  multi-element  commercial 
probes  in  a  transmission  geometry.  The  characterization 
technique  is  based  on  an  acoustic  propagation  model  that 
relates  pair-wise  voltage  measurements  between  multiple 
transducer  elements.  We  use  this  propagation  model  to  de¬ 
rive  a  fast,  semi- analytic  nonlinear  optimization  procedure 
to  estimate  the  transducer  model  parameters  from  self¬ 
measurements.  This  eliminates  the  need  for  external  field 
probes  that  are  different  from  the  sources  themselves,  and 
is  the  first  such  attempt  within  this  formalism  in  acoustics 
or  electromagnetics. 

In  Section  II,  we  derive  the  acoustic  formalism,  which 
is  analogous  to  our  work  in  [13]— [15] .  First,  we  explain  the 
transducer  source  model,  which  is  based  on  the  acous¬ 
tic  source-scattering  matrix  and  derive  the  propagation 
model.  We  further  augment  the  transducer  model  for  the 
expected  electrical  driving  frequencies.  The  transducer 
model  is  then  used  to  modify  the  acoustic  volume  integral 
equations  of  the  inverse  scattering  algorithm.  In  Section 
III,  we  discuss  source  characterization  to  motivate  and 
then  derive  the  optimization  routine  for  the  self-calibra¬ 
tion  technique.  Section  IV  describes  our  experimental  set¬ 
up  and  characterization  results,  including  the  assumptions 
necessary  to  carry  out  the  analysis.  Finally,  in  Section  V, 
the  inverse  scattering  algorithm  is  tested  and  preliminary 
images  of  simple  targets  are  presented  and  discussed. 

Several  provisions  were  required  to  test  this  method 
in  experiment,  which  we  explain  throughout  Sections  IV 
and  V.  We  use  two  linear- array  commercial  hand  probes. 
This  limits  the  setup  to  a  planar  transmission  geometry. 
We  use  a  Verasonics  data  acquisition  system[AU1 :  Please 
provide  manufacturer  name  and  location.],  which  re¬ 
quires  several  assumptions  about  the  nature  of  the  trans¬ 
mit  pulse.  This  leads  to  the  use  of  virtual  sources.  Finally, 
the  implemented  geometry  is  static,  so  test  objects  being 
imaged  must  be  rotated  to  obtain  multiple  transmitter 
views.  We  use  a  robotic  arm  to  rotate  the  objects,  while 
aligning  with  reflection  measurements. 


Then  the  acoustic  propagation  model  is  derived,  which 
is  the  basis  for  the  self-characterization  method.  Finally, 
we  modify  the  volume  integral  equations  so  that  they  are 
consistent  with  the  transducer  model. 


A.  Traditional  Volume  Integral  Equations 


The  acoustic  VIE  for  the  total  pressure  field  in  the 
presence  of  an  inhomogeneous  distribution  of  compress¬ 
ibility,  density,  and  compressive  loss,  is  given  by 


</>(r)  =  Anc(r)  +  ko J g( r,r')  8n{ r')  +  <P(T')dV> 

+  fff(r,r')V'  ■  8p-\v’)V'4>(v')dV' 

k 2  =  k2n 


KnUl 


(1) 

(2) 


where  k20  =  is  the  lossless  background  wave  number 

and  g(r,r')  is  the  free-space  Green’s  function  with  wave 
number  k.  The  incident  field,  0inc(r),  is  the  field  in  the 
absence  of  the  object.  We  define  the  following  contrast 
functions  for  the  density,  compressibility,  and  compressive 
loss,  respectively,  as 


V-'(D  =  $5  -  1 

(3) 

(4) 

Kjq 
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where  k0,  and  a0  are  the  background  material  con¬ 
stants.  The  quantities  5p~l  and  5k  are  unitless,  and  5a 
is  an  absolute  measure  of  compressive  loss  with  units  of 
Pa-1s-1. 

Defining  the  scattered  field  as 

<?Va(r)  =  <Kr)  -  Anc(U  (6) 

we  can  write  the  VIE  for  the  scattered  field  concisely  as 

Aca(r)  =  f  d(r,  r,)<Ar'V(r')d  V',  (7) 


II.  Formulation 

Inverse  scattering  algorithms  are  typically  built  with 
two  volume  integral  equations  (VIEs).  To  use  the  VIEs  in 
experimental  inverse  scattering,  we  must  model  the  fields 
radiated  and  received  by  the  transducers,  and  incorporate 
this  model  into  the  VIEs.  What  follows  is  the  acoustic  an¬ 
alog  of  our  previous  work  [13]— [15] ,  in  which  we  combined 
an  antenna  model  and  electromagnetic  inverse  scattering 
algorithm. 

We  first  give  the  traditional  acoustic  volume  integral 
equations.  We  then  summarize  the  source-scattering  ma¬ 
trix  formulation  and  adapt  it  to  ultrasound  frequencies. 


where  we  define  the  object  function  O(r)  as  the  linear 
operator 


0(r)  =  k20  \  8k( r)  +  j  ]  +  k2V'  ■  8p-l(r)V'.  (8) 


In  the  context  of  inverse  scattering,  (1)  represents  the 
solution  to  the  forward  scattering  problem  and  (7)  is  used 
to  relate  the  material  contrasts  to  scattered  field  mea¬ 
surements  outside  the  object  domain.  Traditionally,  these 
equations  are  used  as  is  to  develop  inverse  scattering  al¬ 
gorithms. 
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B.  Source  Model 


1 )  Field  Expansion:  The  pressure  field  about  an  acous¬ 
tic  source  can  be  expanded  in  terms  of  incoming  and  out¬ 
going  wave  functions  as 

oo  l 

(^/m^/m(r)  T  i  (^) 

l=0m=—l 

where  aim  and  6/m  are  the  multipole  coefficients  for  incom¬ 
ing  and  outgoing  harmonics,  respectively,  r  is  the  position 
vector,  and  K  means  the  regular  part  of  the  corresponding 
Bessel  function.  The  double  sum  is  hereafter  abbreviated 
^2lm  with  the  same  limits.  The  free-space  scalar  wave 
functions  are  given  by  [17] 

^m(r)  =  zi(kr)Ylm(6,4 >),  (10) 

where  Z/(  x)  is  any  one  of  yl  x) ,  h\l\x),  h\2\x),  and 

Yim{0^)  are  the  spherical  angular  harmonics. 

2)  Source  Scattering  Matrix  Formulation :  Following  the 
source  scattering- matrix  formulation  [18],  [19],  a  trans¬ 
ducer  is  connected  through  a  shielded  transmission  line 
to  a  voltage  source  or  receiver,  Fig.  1.  We  define  a0  and 
b0  as  complex  outgoing  and  incoming  modes,  respectively, 
on  the  transmission  line  measured  at  a  reference  plane. 
We  relate  these  modes  to  the  multipole  coefficients  of  the 
pressure  field  through  a  linear  model 

b0  =  IX  +  ywima/m  (11) 

lm 

b lm  ^Irrfio  T  ^  ^ IrriJ'm'^ I'm' (1^) 
I'm' 


where  tim  and  uim  are  transmit  and  receive  coefficients,  re¬ 
spectively.  The  matrix  Sim  j>m>  captures  the  passive  scatter¬ 
ing  properties  of  the  acoustic  source  and  Y  is  the  reflection 
coefficient  seen  looking  into  the  device. 

From  electro-acoustic  reciprocity,  the  following  rela¬ 
tions  can  be  derived  between  the  transmit  and  receive 
coefficients  [18], 


Pressure  Refd 
expansion 
Coefficients 


Fig.  1.  Acoustic  source  and  transmission  line  setup.  a0  and  b0  are  outgo¬ 
ing  and  incoming  modes  on  the  transmission  line,  respectively. 


nlm  =  A(-!)V  (13) 

Zp0CK 

=  (“I  r+mA_ra,iW,  (14) 

where  p0,  c,  and  k  are  the  background  density,  sound 
speed,  and  wave  number,  respectively,  and  Z0  is  the  char¬ 
acteristic  impedance  of  the  line. 

3)  Relation  to  Total  Line  Voltage:  The  use  of  left  and 
right  traveling  line  voltages  is  advantageous  for  microwave 
problems,  because  it  allows  the  source-scattering  formula¬ 
tion  to  be  naturally  adapted  to  microwave  S-parameters 
[13].  Given  the  long  electrical  wavelengths  at  ultrasound 
operating  frequencies  (1  to  10  MHz),  full  network  anal¬ 
ysis  is  not  usually  required.  We  measure  only  the  total 
line  voltage,  so  there  is  no  need  to  define  a  phase  refer¬ 
ence  plane  on  the  transmission  line.  However,  the  electri¬ 
cal  impedance  (and  thus  the  reflection  coefficient)  of  the 
acoustic  source  does  affect  the  system  response.  Here,  we 
modify  the  source  scattering  formulation  for  use  with  the 
total  line  voltage. 

The  total  line  voltage  is  given  as  V0  =  a0  +  b0.  If  the 
acoustic  device  is  purely  receiving,  we  expect  no  input 
mode,  a0  =  0,  so  b0  =  V0 ,  which  we  define  as  the  total 
output  voltage  on  the  line,  Fout: 

Vout  K  'y  'u lrna im  •  (1^) 

lm 


If  the  device  is  purely  transmitting  and  no  other  device 
is  transmitting,  so  there  are  no  other  incoming  waves  in 
the  frame  of  the  transmitter  (ajm  =  0),  then  b0  =  Ya0  and 
bim  =  timao •  Defining  the  total  line  voltage  now  as  Fin,  we 
have 


blm  t lm  ^  _|_  p  (15) 

It  is  enough  to  lump  the  effects  of  the  reflection  coef¬ 
ficient  into  the  reciprocity  relations.  The  source  model  is 
now 


Lout  : 

-  TWin 

(17) 

lm 

blm 

tlmY^im 

(18) 

^ lm 

cr(— l)”A,-m 

(19) 

cr  = 

Z<k  1  +  r) 

2  P0ck2 

(20) 

where  we  define  cr  as  the  reciprocity  constant. 

The  overarching  assumption  is  that  the  devices  do  not 
transmit  and  receive  simultaneously,  or,  if  they  do,  the 
signals  can  be  separated,  for  example,  by  time  gating. 
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It  may  be  difficult  or  impossible  to  measure  Z01  T,  and 
possibly  the  other  coefficients  accurately  (especially  with 
commercial  ultrasound  probes).  However,  we  see  that  u\m 
and  bm  are  simply  linearly  proportional  by  systematic 
constants.  The  important  part  of  the  reciprocity  relation 
is  the  dependence  between  uim  and  tim  on  —m. 

C.  Acoustic  Propagation  Model 

Next  we  derive  an  acoustic  propagation  model  relating 
the  multipole  fields  of  two  transducers  to  voltage  measure¬ 
ments  between  them.  This  model  will  be  used  later  for 
the  characterization  procedure.  An  analogous  propagation 
model  for  antenna  electric  fields  and  microwave  S-param- 
eter  measurements  can  be  found  in  [13]. 

Let  there  be  two  acoustic  transducers,  which  could  be 
elements  on  different  scan  heads,  one  transmitting  in  a 
frame  z,  the  other  receiving  in  a  frame  j.  Both  have  their 
own  field  expansion  and  set  of  transmit  and  receive  co¬ 
efficients.  Assuming  no  multiple  scattering  between  the 
sources,  the  fields  in  each  frame  are  expanded  as 

<KTi)  =  (21) 

lm 

0(rj)  =  (22) 

lm 

where  r,  and  r,  are  frame-centered  position  vectors.  From 
the  addition  theorem,  the  multipole  coefficients  of  the  two 
frames  are  related  by[AU2:  Please  check  all  equation 
numbers  and  equation  citations  carefully.  Manuscript 
had  a  blank  line  numbered  (23).] 

ai'm'  =  TW'm'itAm!  (23) 

lm 

where  is  the  translation  matrix  [17],  and  is  evalu¬ 

ated  with  the  vector  that  points  from  frame  i  to  frame  j. 
Using  (17),  (18),  and  (23),  the  input/output  voltages  of 
the  two  transducers  are  related  by 

VL  =  im-  (24) 

I'm'  lm 

Substituting  (19[AU3:  into  (24)?])  and  taking  m!  — > 
—  m\  we  have  the  acoustic  propagation  model  in  terms  of 
only  transmit  coefficients: 

VL  =  (25) 

I'm'  lm 

where  depends  on  the  properties  of  the  receiver.  Eq. 
(25)  relates  the  transmit  coefficients  of  two  transducers  to 
voltage  measurements  between  them. 


D.  Modified  Acoustic  Volume  Integral  Equations 

Before  using  the  VIEs  in  the  inverse  scattering  algo¬ 
rithm,  we  must  make  them  consistent  with  the  transducer 
model.  Like  [14],  in  which  we  derived  a  new  kernel  to 
directly  link  the  electrical  properties  of  an  object  to  micro- 
wave  S-parameter  measurements,  we  can  derive  an  analo¬ 
gous  kernel  for  the  acoustic  problem.  This  kernel  directly 
links  the  object  acoustic  properties  to  measured  voltages. 
The  kernel  will  be  a  single- argument  scalar  function  spe¬ 
cific  to  the  receiver,  which  will  replace  the  two-argument 
scalar  Green’s  function.  By  reciprocity,  we  will  see  that 
it  is  proportional  to  the  transducer  incident  field,  imply¬ 
ing  that  only  the  incident  field  is  required  to  calibrate  an 
acoustic  inverse  scattering  experiment. 

1 )  Normalized  Fields:  We  first  define  the  normalized  in¬ 
cident  pressure  field  by  substituting  the  transducer  model 
into  (9)  when  the  transducer  is  purely  transmitting. 

Ur)  =  (26) 

lm 

=  Un  (27) 

lm 

=  Vjjfi,  (28) 

where  0inc(r)  is  the  normalized  incident  field.  It  is  the  field 
given  by  only  the  transmit  coefficients. 

We  define  the  normalized  total  field, 

4>(r)  =  <p(r)/ Vin.  (29) 

This  is  the  same  total  field  given  by  (1)  but  due  to  the 
normalized  incident  field. 

2)  Receiver  Kernel:  To  derive  the  receiver  kernel,  we 
first  expand  the  scattered  field  from  the  VIE  as  incoming 
waves  in  the  frame  of  the  receiver.  After  using  the  trans¬ 
ducer  model,  we  rearrange  the  expression  to  resemble  a 
VIE  from  which  we  can  identify  the  new  kernel. 

The  addition  theorem  for  the  scalar  Green’s  function  is 
given  by  [17], 

g(  r,r7)  =  Ay}R^Jr)V4(r'),  (30) 

lm 

where  we  have  taken  the  case  |r|  <  |r'|.  The  wave  function 
fj*  is  the  same  as  ?/>,  but  with  the  angular  harmonics  con¬ 
jugated.  Substituting  (30)  into  (7)  and  rearranging,  we 
can  write  the  scattered  field  in  the  frame  of  a  receiver  in 
terms  of  incoming  waves: 

<?Va(r)  =  (3U 

lm 


where  the  multipole  coefficients  are 
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aim  =  ikJ (32) 

Substituting  (32)  into  (17),  we  can  write  the  output 
voltage  as 

Uut  =  v’yUv'AV' .  (33) 

lm 

By  rearranging  this,  we  can  identify  the  receiver  kernel, 
which  we  label  g(vr), 

Lout  =  f 9(r')0( r')#-')d  V'  (34) 

g(v')  =  ikjyimiptm(T').  (35) 

lm 

The  function  g(rf)  directly  relates  the  output  voltage  of 
the  receiver  to  the  material  properties  of  the  object 
through  the  transducer  model.  It  is  the  acoustic  analog  of 
the  receiver  kernel  derived  in  [14]  for  electromagnetics. 

Finally,  let  there  be  two  transducers,  one  transmitting 
in  frame  z,  the  other  receiving  in  frame  j.  The  total  pres¬ 
sure  field  in  the  object  is  due  to  the  transmitter.  Dividing 
both  sides  of  (34)  by  the  transmit  voltage,  we  can  write 
the  scattered  field  VIE  in  terms  of  the  normalized  total 
field  of  the  transmitter  and  the  input / output  line  voltages 
as 

//  =-  V(r%(r')dV',  (36) 

where  g  j(r')  is  the  kernel  for  the  receiver. 

3)  Reciprocity  Relation:  By  reciprocity,  we  can  show 
that  g r')  is  related  to  the  incident  field  of  the  receiver. 
Substituting  the  transmit  and  receive  coefficient  reciproc¬ 
ity  relations,  (19),  into  (35), 

g{v’)  =  (37) 

lm 

and  using  the  relation  Y*m  =  (— 1  )mYi_m  in  the  wave  func¬ 
tion,  we  get 

g( r')  =  ikcTY)i,m^im(r '),  (38) 

lm 

which  is  the  expansion  for  the  normalized  incident  pres¬ 
sure  field  multiplied  by  a  scaling  factor 

g(r')  =  ikc^iDC( r').  (39) 

Finally,  we  write  the  scattered  field  transmit /receive 
voltage  ratio  in  terms  of  the  normalized  incident  field  of 
the  receiver  and  normalized  total  field  of  the  transmit¬ 
ter  as 


1 m 

This  equation  consistently  links  the  material  properties 
we  wish  to  image  to  the  voltages  we  measure.  Only  the 
normalized  incident  fields  of  the  transmitters  and  receiv¬ 
ers  are  required  to  do  this.  We  will  use  (40)  in  the  inverse 
scattering  algorithm  when  comparing  forward  model  pre¬ 
dictions  to  measurements.  This  is  the  acoustic  analog  of 
the  expression  derived  in  [14]  for  microwave  S-parameter 
measurements. 


III.  Ultrasound  Source  Characterization 

The  next  step  is  to  determine  the  incident  fields  of  the 
transducers,  which  are  the  quantities  needed  to  character¬ 
ize  an  inverse  scattering  experiment.  The  incident  fields 
can  be  measured  directly  or  computed  with  (27)  after  the 
transmit  coefficients  are  found.  The  transmit  coefficients 
can  be  estimated  with  simulation  or  determined  from 
measurement. 

We  would  like  to  determine  the  transmit  coefficients 
in  a  way  that  stays  within  the  formalism  of  the  source¬ 
scattering  model  because  of  the  rigor  and  versatility  of  the 
preceding  derivation.  Here,  we  discuss  several  aspects  of 
this  problem  and  the  motivation  for  our  method  of  choice 
given  the  experimental  circumstances.  Keep  in  mind  that 
we  intend  to  characterize  the  transducer  elements  of  com¬ 
mercial  hand  probes. 

A.  Hydrophone  Characterization 

If  the  pressure  field  from  an  acoustic  source  can  be 
sampled  with  a  calibrated  hydrophone,  then  it  is  possible 
to  estimate  the  transmit  coefficients  by  solving  a  linear 
inverse  problem.  The  source  model  is  given  by 

0(**n)  ^jJlm^Plm^n)i  (^1) 

lm 

where  rn  are  the  locations  of  N  hydrophone  pressure  field 
measurements,  and  Vin  is  the  frequency  component  of  the 
driving  voltage.  With  enough  measurements  for  a  given 
frequency,  tim  is  estimated  by  solving  a  linear  system  of 
equations  over  (41). 

Although  direct  pressure  measurements  lead  to  the 
simplest  estimation  of  the  transmit  coefficients,  the  main 
challenges  with  this  approach  are  1)  alignment  and  posi¬ 
tioning,  2)  obtaining  the  hydrophone  transfer  function, 
and  3)  repeatability.  Furthermore,  the  pattern  and  fre¬ 
quency  response  of  the  hydrophone  must  be  taken  into 
account. 

We  are  interested  in  characterizing  multi-element  com¬ 
mercial  transducer  probes.  Commercial  probes  present  a 
challenge  to  the  method  described  here  because  of  the  dif- 
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ficulty  of  positioning  individual  elements  as  well  as  model¬ 
ing  the  system  electronics. 

B.  Simulation 

In  our  previous  work  characterizing  antennas  [13],  we 
took  advantage  of  the  fact  that  we  could  make  our  an¬ 
tenna  construction  and  its  CAD  model  for  our  commercial 
simulation  package  (Ansoft  HFSS)  completely  consistent. 
This  enabled  us  to  estimate  the  antenna  transmit  coef¬ 
ficients  using  simulated  electric  fields  and  an  expression 
similar  to  (41).  However,  complications  arise  when  trying 
to  simulate  ultrasound  transducers  for  the  same  purpose. 

First,  we  do  not  fabricate  ultrasound  transducers  our¬ 
selves.  They  are  an  assortment  of  piezoelectric  resonators, 
matching  layer  (to  water),  impedance  matching  (to  driv¬ 
ing  circuit),  grounding  layer,  plastic  membrane,  water 
proofing,  and  focusing  lens  [20].  It  is  possible  to  obtain 
a  rough  mechanical  layout  by  obtaining  an  X-ray  image 
of  a  transducer,  but  the  exact  material  properties  and 
boundary  conditions  are  difficult  to  obtain.  Even  if  all 
of  the  mechanical  properties  are  known,  simulating  the 
electro-acoustic  response  requires  a  multi-physics  simula¬ 
tion  package  which  can  model  anisotropic  piezoelectric 
material. 

Second,  the  open-access  package  Field  II  [21]  simulates 
the  free-space  radiation  patterns  of  transducer  apertures 
and  arrays.  However,  it  is  based  simply  on  the  linear  im¬ 
pulse  response  of  ideal  apertures,  and  does  not  account 
for  the  driving  circuits  or  any  mechanical  features  of  the 
transducers.  Thus,  the  simulated  fields  are  not  suitable  for 
full-wave  characterization. 

Although  piezoelectric  transducers  are  difficult  to  fully 
simulate,  the  simplicity  of  CMUT  transducers  has  allowed 
successful  electro-mechanical  simulations  which  agree  with 
measurements  [22].  The  transmit  coefficients  of  these  de¬ 
vices  might  be  estimated  from  simulation,  which  we  note 
as  possible  future  work. 

C.  Nonlinear  Inversion  of  the  Propagation  Model 

Without  full  electromechanical  simulation  or  hydro¬ 
phone  measurements  of  our  commercial  transducers,  the 
final  option  is  to  obtain  the  transmit  coefficients  by  in¬ 
verting  the  propagation  model.  This  allows  us  to  stay 
within  the  formalism  of  both  the  source-scattering  model 
and  modified  VIEs,  and  estimate  the  transmit  coefficients 
directly  from  measurement.  It  represents  a  hydrophone¬ 
less,  multi-source  self-calibration  technique.  Two  or  more 
transducers  are  required,  with  known  locations  and  mea¬ 
sured  transmit  and  receive  line  voltages,  from  which  pair¬ 
wise  transducer  responses  are  used  to  build  an  inverse 
problem  around  (25). 

The  advantage  of  this  approach  is  that  all  the  necessary 
physics  are  contained  in  the  propagation  model:  transla¬ 
tion  matrices  capture  wave  propagation,  and  reciprocity 
relates  the  transmit  coefficients  of  different  transducers 
to  each  other.  The  need  for  multiple  pair-wise  measure¬ 


ments  is  satisfied  by  imaging  setups  with  multi-element 
commercial  probes.  The  disadvantage  is  that,  because  the 
transmit  coefficients  inherently  capture  the  absolute  phase 
of  the  radiated  fields,  the  locations  of  the  transducer  refer¬ 
ence  frames  must  be  known  or  assumed  with  submillime¬ 
ter  accuracy.  We  must  also  sample  the  fields  throughout 
the  region  in  which  we  intend  to  compute  the  incident 
field. 

The  propagation  model  is  a  nonlinear  function  of  the 
transmit  coefficients,  but  only  polynomial  to  second  order. 
The  weakly  nonlinear  form  allows  us  to  derive  a  semi- 
analytic  conjugate  gradient  search  for  the  minimum  of  a 
least-squares  function.  This  is  advantageous  for  large  data 
sets  and  fast  optimizations.  We  develop  this  in  the  next 
section  and  will  use  it  to  estimate  the  transmit  coefficients 
from  measurements. 

1 )  Cost  Function  and  Forward  Model:  We  estimate  the 
transmit  coefficients  by  minimizing  a  cost  function  com¬ 
paring  measured  voltages  to  propagation  model  predic¬ 
tions. 

Following  [23] ,  the  least  squares  cost  function  for  a  gen¬ 
eral  nonlinear  problem  is 

2F(m)  =  ||  g(m)  -  d||*  +  ||m  -  (42) 

The  vector  norms  are  defined  over  model  and  data 
spaces  through  their  respective  inverse  covariance  matri¬ 
ces,  C^1  and  Cfi1.  The  elements  of  the  forward  model  vec¬ 
tor,  g(m),  and  data  vector,  d,  are  given  by 

gjX m)  =  (43 ) 

I'm'  Im 

vK 

d»  =  IM,  (44) 

^in 

with  the  self  term  excluded.  The  vector  of  model  param¬ 
eters  m  contains  the  transmit  coefficients  of  all  N  trans¬ 
ducers,  written  in  vector  form  over  all  harmonics  as 

(45) 


In  the  case  in  which  all  the  transducers  in  the  problem 
are  identical,  the  model  has  quadratic  form,  and  we  es¬ 
sentially  solve  the  problem  c  =  x2.  However,  for  multiple 
distinct  transducer  elements,  for  which  each  transducer 
has  its  own  set  of  transmit  coefficients,  the  forward  model 
takes  the  form  c  =  xy.  In  the  former  case,  the  answer  is 
unique  to  within  a  sign.  For  the  latter,  it  is  impossible 
to  uniquely  separate  the  product  of  two  unconstrained 
variables. 

The  properties  of  the  forward  model  suggest:  1)  that 
we  should  attempt  to  solve  a  self-characterization  problem 
only  when  we  can  assume  all  the  transducers  are  identi- 
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cal,  and  2)  the  polynomial  form  enables  us  to  derive  an 
analytical  expression  for  the  step  length  in  a  conjugate 
gradient  minimization,  enabling  a  fast  gradient  search. 

To  facilitate  the  derivations  here  and  in  the  appendix, 
we  write  the  forward  model  in  matrix  notation  as 

gji( m)  =  tjNajjtj,  (46) 

where  ()*  is  vector  transpose,  N  is  a  diagonal  matrix 
containing  cr(— l)m,  and  the  rows  of  are  evaluated  at 
m!  — ►  —m'. 


2)  Gradients:  Assuming  the  data  are  independent,  the 
gradient  of  the  cost  function  with  respect  to  the  model 
parameters  is  given  by 


=  E 


df_ 

<9m 


-^(gjXm)-dji)  +  C  (47) 

at 


where  dg^/dm  is  the  vector  of  partial  derivatives  of  the 
forward  model  with  respect  to  the  model  parameters. 

The  partial  derivatives  of  (47)  are  computed  with  the 
help  of  these  matrix  identities: [AU4:  Is  the  superscript 
capital  T  supposed  to  be  a  superscript  lowercase  t  to 
indicate  a  vector  transpose  as  in  (46)?] 


<9xtBx 

dx 

<9xTa 

dx 


--  (B  +  Bt)x 
<9a  Tx 

“&T  =  a‘ 


(48) 

(49) 


(^M  {if n  Tn— l)?7n) 

(^M  T  n— V>  T  n— l) 


(55) 


where  (,)  is  a  simple  dot  product. 

The  step  length  an  is  found  by  minimizing  the  cost 
function  at  each  step.  Given  the  form  of  the  forward  mod¬ 
el,  this  can  be  done  analytically  as  described  in  the  ap¬ 
pendix. 

4)  Summary  of  Optimization  Routine:  We  derived  a 
nonlinear  optimization  routine  to  invert  the  acoustic  prop¬ 
agation  model  for  the  transducer  transmit  coefficients  giv¬ 
en  pair-wise  voltage  responses.  All  of  the  preceding  steps 
were  used  to  implement  this  routine.  Experimental  results 
follow  in  the  next  section. 

Much  time  was  spent  analyzing  this  inversion  routine  in 
simulation.  Results  obtained  in  our  simulations  were  con¬ 
sistent  with  the  algorithm  presented  here.  With  perfect 
data  and  perfect  knowledge  of  the  source  locations,  the 
transmit  coefficients  for  multiple  identical  sources  could 
be  recovered  uniquely  to  within  a  sign,  from  a  wide  range 
of  initial  conditions.  It  was  also  determined  that  when 
solving  for  multiple  different  transducer  transmit  coeffi¬ 
cients,  that  the  synthetic  measured  data  could  be  fit  quite 
well,  but  the  solution  was  not  correct  and  dependent  on 
the  initial  condition,  demonstrating  non-uniqueness. 


IV.  Ultrasound  Probe  Characterization 
Experiment 


If  the  transducers  are  identical,  t  =  t j 
and  j,  then 

=  U,  for  all  i 

Baji 

-Jj-  =  (Na  ji  +  (Na^t. 

(50) 

If  the  transducers  are  different,  then  the  gradients  are 

a 

& 

II 

(51) 

d9ji  t»t  , 

<9t.  “  NaA- 

(52) 

3)  Conjugate  Gradient  Updates:  Regardless  of  whether 
we  are  solving  for  identical  or  different  transducer  trans¬ 
mit  coefficients,  the  conjugate  gradient  updates  are  de¬ 
fined 


=  mn-l  -  (53) 

vn  =  7  n  +  PnVn-l,  (54) 

where  an  is  the  step  length  and  is  the  steepest  decent 
vector.  j3n  is  given  by 


To  test  the  nonlinear  source  characterization  in  the  pre¬ 
vious  section,  we  built  the  setup  shown  in  Fig.  2.  This  setup 
was  also  used  for  the  inverse  scattering  (see  next  section) 
so  that  the  characterization  and  imaging  were  performed 
in  identical  geometries.  Two  commercial  Phillips  ATL 
L7-4  transducers[AU5:  Please  provide  manufacturer 
name  and  location.]  were  mounted  facing  each  other  in 
a  water  tank.  Each  probe  had  128  transducer  elements  in 
a  linear  array.  The  probes  were  connected  to  a  Verason- 
ics  data  acquisition  system,  which  was  programmed  to 
collect  all  combinations  of  transmit  and  receive  signals 
between  the  opposite-side  elements.  The  mount  was  cus¬ 
tom  made  out  of  king  starboard. [AU6:  Please  provide 
manufacturer  name  and  location.]  The  mounting  holes 
are  separated  by  1  cm  to  allow  linear  translation  of  the 
probes  relative  to  one  another.  The  blue  holders  were  cre¬ 
ated  with  the  help  of  the  University  of  Michigan  3-D  Lab. 
The  probes  were  laser  scanned  and  the  scan  was  used  to 
create  a  3-D  CAD  model  of  the  holders.  The  holders  were 
then  fabricated  with  a  3-D  printing  system. 

We  took  transmission  measurements  in  probe  positions 
that  would  sweep  out  the  imaging  domain.  This  allowed  us 
to  sample  the  field  both  at  points  in  the  imaging  domain 
as  well  as  at  the  receiver  locations.  This  setup,  however, 
restricts  our  field  samples  to  a  plane,  which  is  a  compro¬ 
mise  we  made  to  avoid  more  involved  3-D  geometries.  The 


IEEE  TRANSACTIONS  ON  ULTRASONICS,  FERROELECTRICS,  AND  FREQUENCY  CONTROL,  VOL.  TBC,  NO.  TBC,  TBC  TBC 


Fig.  2.  Photo  of  setup  for  nonlinear  source  characterization  and  inverse 
scattering  experiment.  A 


transducer  elements  of  the  L7-4  probe  are  known  to  be 

4  mm  tall,  0.3  mm  wide,  with  a  spacing  of  0.4  mm. 

Transmitted  signals  were  two  cycle  pulses  with  a 

5  MHz  center  frequency,  which  was  the  center  frequency 
used  for  the  position  optimization  in  the  next  section.  The 
Verasonics  system  was  set  to  sample  at  4  samples  per 
wavelength,  giving  a  sample  rate  of  20  MHz.  In  water,  the 
wavelength  at  5  MHz  is  300  pm.  The  characterization  and 
imaging  were  performed  at  3.75  MHz,  where  the  wave¬ 
length  is  400  pm  in  water. 

A.  Position  Inversion  and  Virtual  Sources 

Before  performing  the  characterization,  we  must  know 
or  assume  the  locations  and  orientations  of  the  transducer 
reference  frames.  This  is  not  easily  done  mechanically  be¬ 
cause  of  the  probably  imprecise  positioning  of  the  trans¬ 
ducer  elements  in  the  transducer  housing,  so  we  use  time 
of  flight.  The  time  of  flight  is  recorded  between  all  (or 
many)  elements  on  the  opposing  probes  and  used  to  solve 
a  triangulation  problem  for  the  position  and  orientation  of 
one  array  relative  to  the  other.  We  assume  we  know  the 
spacing  of  the  array  elements  and  that  all  the  elements  lie 
on  a  line. 

We  mention  several  experimental  conditions.  We  use 
time  zero  of  the  time  traces  as  a  common  reference,  be¬ 
cause  the  receive  channels  of  the  Verasonics  begin  record¬ 
ing  at  the  same  time  on  each  occurrence  of  a  transmit 
pulse.  The  identical  timing  of  the  channels  is  the  only 
common  reference  between  all  transmit  and  receive  data. 
Second,  this  allows  measurements  of  the  delay  in  firing  of 
the  transmit  pulse  equivalent  to  about  one  wavelength. 
Third,  each  transducer  has  a  lens  on  its  face  with  a  lower 
speed  of  sound.  This  lens  is  designed  to  focus  the  radiation 
from  each  element  in  the  elevation  direction,  so  that  fo¬ 
cused  beams  used  for  everyday  diagnostic  ultrasound  are 
restricted  to  a  thin  plane  out  to  a  modest  distance  from 
the  transducer.  This  lens  introduces  a  delay  on  both  the 
transmit  and  receive  signals.  Last,  the  pulse  may  be  de- 


Fig.  3.  Envelope  of  time  trace.  The  time  of  flight  is  determined  by  the 
time  of  the  envelope  maximum  less  half  the  pulse  width, 

layed  either  by  the  system,  the  cables,  or  probe  head,  and 
we  may  not  be  able  to  quantify  this  delay. 

We  measured  the  time  of  flight  relative  to  the  begin¬ 
ning  of  the  time  trace  using  an  envelope-max,  shown  in 
Fig.  3.  The  firing  delay  was  accounted  for  by  subtracting 
half  a  pulse  length  in  water  from  the  time  of  flight  mea¬ 
sured  at  the  peak  of  the  envelope.  To  properly  account  for 
the  delay  induced  by  the  lenses,  one  approach  is  to  solve 
Fermat’s  least  action  integral  for  the  two-layered  bound¬ 
ary,  which  we  did  not  do  here  because  the  properties  of 
the  lens  were  not  known. 

The  effects  of  not  accounting  for  the  lenses,  or  other 
delays,  in  a  time  of  flight  position  estimation  are  that 
we  will  find  the  location  of  a  virtual  source  behind  the 
actual  source  location.  We  rationalize  this  as  follows:  1) 
accounting  for  the  lens  in  the  propagation  model  inversion 
is  not  trivial,  and  2)  if  the  lenses  do  not  adversely  affect 
the  radiation  pattern  in  the  azimuth  plane  (the  plane  of 
the  array),  then,  using  the  virtual  positions,  we  can  solve 
for  the  transmit  coefficients  of  equivalent  sources  at  the 
virtual  locations.  In  doing  so,  the  effects  of  the  lenses  are 
accounted  for  by  a  free-space  spatial  offset.  Finally,  if  we 
are  finding  virtual  source  locations  and  transmit  coeffi¬ 
cients,  we  can  now  also  lump  the  effects  of  miscellaneous 
system  delays  into  the  virtual  source  location.  This  was 
the  reasoning  going  forward  and  the  assumptions  under 
which  we  performed  the  following  computations  and  ex¬ 
periments.  Fig.  4  shows  a  diagram  of  the  transition  from 
what  we  believe  to  be  the  actual  system  to  the  virtual 
source  configuration. 

The  forward  model  for  the  position  inversion  is  dia¬ 


grammed  in  Fig.  5  and  is  given  by 

hi  =  rijl c  (56) 

n,-  =  |-  (Cl  +  d*Vi)  +  C2  +  djv2\  (57) 

dt  =  ((128  +  l)/2  -  i)d,  i  =  1...128  (58) 

dj  =  ((128  +  l)/2  -  j)d,  j  =  1. . .128,  (59) 


where  tLJ  is  the  time  of  flight  between  two  elements,  c  is 
the  speed  of  sound  in  water,  and  is  the  distance  be¬ 
tween  two  elements.  The  vectors  ci  and  C2  point  to  the 
centers  of  each  array.  The  vectors  di\l  and  d-v2  point  to 
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Assumed  virtual  source  and  receiver 
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Fig.  4.  Actual  and  virtual  sources  and  receivers.  The  signal  fires  at 
emanates  from  the  transducer  at  q,  arrives  at  the  receiving  transducer  at 
^2,  and  is  recorded  at  t3.  The  virtual  source  and  receiver  lump  any  delays 
into  a  spatial  offset.  A. 


the  locations  of  the  array  elements,  with  d  being  the  sepa¬ 
ration  of  the  elements. 

B.  Speed  of  Sound  in  Water 

The  speed  of  sound  of  the  water  bath  was  determined 
by  the  temperature  at  the  time  the  data  were  taken 
and  the  lookup  table  in  [24].  This  was  periodically  cross 
checked  with  pulse-echo  reflections  off  a  plate  for  a  trans¬ 
ducer  stepped  on  a  micrometer.  The  two  methods  agreed 
to  within  0.1%  at  3.75  MHz  for  speeds  around  1495  m/s 
at  24.3°C. 

C.  Position  Inversion  Results 

After  determining  the  times  of  flight  between  all  ele¬ 
ments,  we  minimized  an  LI  cost  functional  over  (56);  the 
LI  norm  was  used  to  automatically  eliminate  outliers  in 
the  time  of  flight  data  created  by  the  envelope  search. 

We  solved  for  the  location  and  orientation  of  the  second 
transducer  relative  to  the  first.  Thus,  we  let  cq  =  0  and 
v1  =  x  and  the  six  model  parameters  are  then  the  ele¬ 
ments  of  the  vectors  c2  and  v2.  The  initial  conditions  used 
were  c2?0  =  c  *  £centery,  where  £center  is  the  time  of  flight 
between  the  two  center  elements  and  v2  =  — x,  which  as¬ 
sumes  that  the  arrays  are  initially  parallel  and  centered. 

This  optimization  was  performed  for  seven  positions. 
The  results  are  shown  in  Fig.  6.  The  probe  position  sepa¬ 
rations  were  consistent  with  our  1  cm  spacing,  but  the  re- 


Fig.  5.  Diagram  for  position  estimation  problem.  The  vectors  C2  and  ci 
point  to  array  centers,  V]_  and  v2  point  from  the  array  centers  to  element 
locations,  and  d  is  the  element  spacing. 
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Fig.  6.  Split  plot  of  the  probe  element  positions.  Filled  circles  are  the 
element  locations  of  probe  1  (reference) .  Open  circles  are  the  element  lo¬ 
cations  of  probe  2  for  all  seven  positions  based  on  the  position  inversion. 
The  y- axis  is  compressed. 


ceive  array  was  shifted  along  the  array  direction  by  nearly 
a  millimeter.  This  was  confirmed  visually,  and  corresponds 
to  a  translation  equivalent  to  several  array  elements.  We 
ignored  the  z  coordinate  offset  and  z  orientation  compo¬ 
nent,  because  1)  the  position  inversion  is  not  sensitive  to 
that  dimension  unless  the  array  is  grossly  tilted  in  the  x-z 
plane,  and  2)  these  variations  are  small  compared  with 
variation  in  the  element  radiation  patterns  in  the  eleva¬ 
tion  direction. 

D.  Characterization  Results 

Having  estimated  the  locations  of  the  sources,  we  then 
estimated  the  transmit  coefficients  using  the  method  out¬ 
lined  in  Section  III-C.  We  did  this  in  the  frequency  do¬ 
main  at  3.75  MHz,  which  was  then  the  frequency  used 
for  imaging.  Although  more  frequencies  ultimately  help 
the  inverse  scattering  algorithm,  a  single  frequency  is  suf¬ 
ficient  to  confirm  that  the  characterization  and  imaging 
algorithm  work  together.  This  frequency  was  the  lowest 
at  which  we  could  drive  the  transducers  for  clean  signals; 
lower  frequencies  allow  more  freedom  for  phase  and  posi¬ 
tion  errors. 

Although  the  derivation  of  (40)  is  rigorous,  we  men¬ 
tion  the  problem  of  determining  the  systematic  constants 
in  the  context  of  this  experiment.  These  include  Z0,  T, 
and  the  Verasonics  system  electronics  transfer  functions, 
which  would  convert  sampled  signals  to  absolute  voltages. 
We  did  not  measure  these  for  the  following  reasons.  By 
determining  the  transmit  coefficients  in  the  same  system 
configuration  used  to  take  data  for  the  inverse  scattering, 
the  constants  carry  through  to  the  computation  of  the 
normalized  incident  field  in  the  VIE.  In  other  words,  we 
estimate  scaled  transmit  coefficients  using  digitized  volt¬ 
age  data.  The  same  transmit  coefficients  are  used  to  com¬ 
pute  the  normalized  incident  fields  in  the  VIEs.  To  port 
the  transmit  coefficients  to  another  system,  new  constants 
would  have  to  be  determined. 

Because  we  cannot  take  data  out  of  plane,  we  only  need 
to  describe  the  element  radiation  pattern  in  the  imaging 
plane,  and  so  we  chose  to  invert  for  multipole  harmon¬ 
ics  up  to  l  —  8.  We  assumed  all  transducer  elements  are 
identical. 
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Fig.  7.  Normalized  residual  as  a  function  of  conjugate  gradient  iteration 
number. 

The  frames  of  elements  on  the  second  probe  were  ro¬ 
tated  by  means  of  rotation  matrices  in  the  propagation 
model  [13].  After  many  trials,  we  found  the  best  results  by 
using  data  from  a  subset  of  elements  from  each  probe  and 
a  few  positions.  For  the  presented  examples,  we  used  data 
from  every  tenth  element  and  two  positions  to  obtain  the 
results  below.  These  were  chosen  to  subsample  each  array 
and  avoid  underfitting. 

Fig.  7  shows  the  residual  with  each  iteration  of  con¬ 
jugate  gradient  minimization.  The  residual  drops  mono- 
tonically,  meaning  that  our  algorithm  was  coded  correctly. 
The  behavior  of  the  residual  is  more  complicated  than  for 
a  typical  linear  problem;  the  residual  of  a  linear  problem 
would  look  smoothly  logarithmic. 

Fig.  8  shows  the  real  and  imaginary  parts  of  the  esti¬ 
mated  transmit  coefficients.  This  plot  is  essentially  the 
spectral  power  of  each  multipole  in  the  radiation  pattern 
of  the  elements. 

Fig.  9  shows  the  measured  and  predicted  magnitude 
and  phase  of  the  data  for  a  small  subset  of  transmitter- 
receiver  pairs  after  the  transmit  coefficients  were  found. 
The  groups  are  receivers  per  transmitter.  We  see  that  the 
overall  shape  of  the  magnitude  is  well  fitted,  and  that  the 
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Fig.  8.  Normalized  transmit  coefficients:  (a)  real  part  and  (b)  imaginary 
part.  Horizontal  axis  is  the  (/,  m)  harmonic  index. 


(M 

Fig.  9.  (a)  Magnitude  and  (b)  phase  of  measured  and  predicted  data  for 
a  sampling  of  transmit  and  receive  elements.  The  apparent  groups  in  the 
magnitude  plot  are  receivers  for  a  given  transmit  element. 

phase  is  fit  well,  to  within  20°.  This  indicates  that  the  in¬ 
version  was  working  to  correctly  fit  the  data. 

The  normalized  pressure  field  computed  with  the  trans¬ 
mit  coefficients  is  shown  in  Fig.  10.  This  is  the  normalized 
magnitude  in  the  plane  of  the  probes,  and  is  the  incident 
field  we  will  use  in  the  inverse  scattering  algorithm.  The 
asymmetric  side  lobes  come  from  having  incomplete  infor¬ 
mation  in  those  directions;  the  fields  are  only  constrained 
where  we  took  data  with  the  second  probe. 

V.  Ultrasound  Inverse  Scattering  Experiment 

We  formed  images  of  compressibility  and  absorption 
with  the  acoustic  inverse  scattering  algorithm  in  [9].  The 
algorithm  is  based  on  the  Born  iterative  method  (BIM) 
with  Neumann  series  forward  scattering  solver.  We  re¬ 
placed  the  scattered  field  VIE  in  that  algorithm  with  the 
modified  VIE  of  (40),  and  used  normalized  incident  fields 
computed  with  the  estimated  transmit  coefficients. 

We  formed  images  of  the  2-D  cross  section  of  cylindri¬ 
cal  targets  using  the  full  3-D  inverse  scattering  algorithm 
using  transmission  data  only.  Because  of  the  limited  angle 
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Fig.  10.  Normalized  magnitude  of  acoustic  field  in  the  transducer  plane 
computed  with  the  solved  transmit  coefficients.  Colorscale  is  normlized 
acoustic  field  magnitude  (unitless).  A, 
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Fig.  11.  Photo  of  6-axis  positioning  arm  used  to  place  and  rotate  objects 
and  close-up  of  metal  rod. 

acquisition  geometry  of  the  opposing  probes,  we  rotated 
the  objects  to  adequately  sample  the  scattered  fields.  The 
objects  hung  from  the  metal  rod  and  are  weighed  down. 
We  assumed  the  imaging  domain  was  very  thin  and  ig¬ 
nored  out-of-plane  scattering.  The  modified  algorithm  was 
tested  in  simulation  with  this  geometric  configuration  and 
the  transmit  coefficients  estimated  previously,  which  con¬ 
firmed  that  the  inversion  algorithm  and  forward  solver 
were  self-consistent,  the  details  of  which  are  not  included 
here. 

Several  objects  were  used  to  test  the  inverse  scattering 
algorithm.  The  objects  consisted  of  metal,  plastic,  and 
worm  rubber  cylinders.  The  hard  objects  tested  resolu¬ 
tion.  Worm  rubber  has  a  speed  of  sound  less  than  water 
and  attenuation  greater  than  water  and  was  used  to  test 
the  contrast. 

A.  Rotator  Alignment 

We  rotated  the  objects  using  the  six- axis  robotic  arm 
shown  in  Fig.  11.  Given  the  sensitivity  of  the  algorithm  to 
position  errors,  we  aligned  the  center  of  rotation  as  best 
we  could  at  the  coordinate  center  of  the  object  domain.  A 
smooth  metal  cylindrical  rod  was  attached  to  the  center  of 
the  last  arm  axis.  This  rod  was  used  for  attaching  objects 
and  was  aligned  in  the  x-  and  //-directions  using  pulse  echo 
measurements. 

B.  Image  Reconstructions 

We  formed  images  at  3.75  MHz  using  12  or  36  object 
rotations  and  a  subset  of  transmitters  and  receivers  on 
each  probe.  The  imaging  domain  was  about  2  cm  on  a 
side. 

1)  Example  1:  Worm  Rubber  Cylinder:  An  approxi¬ 
mately  4-mm  diameter  cylinder  of  worm  rubber  was  im¬ 
aged  [Fig.  12(a)].  From  mass/ volume  measurements  of 
the  worm  rubber  samples,  their  density  was  1.02  g/cm3. 
The  speed  of  sound  was  estimated  from  differential  time 
of  flight  to  be  1406  m/s.  From  this,  the  compressibility 
was  determined  to  be  4.98e-10  1/N,  which  equates  to  a 
contrast  with  water  of  0.102.  Attenuation  coefficient  was 
measured  separately  at  3.75  MHz  as  3.2  dB/cm/MHz 
[25].  Twelve  rotations  of  30°  each  were  used  for  this  object. 


X  ion)  X  lcm> 

(b)  (e) 

Fig.  12.  (a)  Photo  of  worm  rubber  target,  (b)  Recovered  compressibility 
contrast,  5k,  (colorscale:  unitless),  and  (c)  absorption  contrast,  5a,  (col- 
orscale:  units  of  compressive  loss  Pa-1s-1). 

Figs.  12(b)  and  12(c)  show  images  formed  by  the  in¬ 
verse  scattering  algorithm  of  compressibility  contrast,  5  k, 
and  absorption  contrast,  5a ,  respectively,  as  defined  by  (4) 
and  (5).  The  inverse  scattering  algorithm  was  terminated 
after  12  BIM  steps  with  12  conjugate  gradient  steps  per  it¬ 
eration.  The  object  location,  shape,  and  approximate  size 
are  recovered.  The  value  of  compressibility  is  too  large  by 
a  factor  of  six  and  the  absorption  shows  a  ring  artifact. 
These  are  most  likely  due  to  the  thin  imaging  domain,  in 
which  the  recovered  values  are  trying  to  compensate  the 
predicted  scattered  fields.  This  effect  is  discussed  in  more 
detail  in  the  next  section.  However,  the  signs  of  the  two 
contrasts  are  correct,  which  indicates  that  (40)  was  de¬ 
rived  correctly.  Similar  effects  are  present  in  the  examples 
that  follow.  The  streaks  along  the  outside  of  the  imaging 
domain  are  rotation  artifacts. 

2)  Example  2:  Two  Metal  Filaments:  We  imaged  two 
thin  metal  filaments  to  test  resolution,  even  though  the 
material  hardness  is  not  captured  by  the  physics  of  the 
inverse  scattering  algorithm.  The  filaments  are  shown  in 
Fig.  13(a).  The  rods  were  brass  with  0.5  mm  diameter  and 
separated  by  approximately  2  mm. 

Figs.  13(b)  and  13(c)  show  images  formed  by  the  in¬ 
verse  scattering  algorithm  of  compressibility  contrast, 
5k,  and  absorption  contrast,  5a,  respectively.  The  inverse 
scattering  algorithm  was  terminated  after  5  BIM  steps. 
The  object  location  and  shape  of  the  object  are  recovered. 
The  compressibility  contrast  is  negative  and  absorption  is 
positive.  This  is  the  correct  direction  for  compressibility, 
because  we  expect  the  compressibility  of  brass  to  be  less 
than  water.  The  ringing  artifacts  and  streaks  are  due  to 
the  use  of  only  twelve  rotations. 
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(a) 


Fig.  13.  (a)  Photo  of  two  metal  filaments,  (b)  Recovered  compressibility 
contrast,  5k,  (colorscale:  unitless),  and  (c)  absorption  contrast,  5a,  (col- 
orscale:  units  of  compressive  loss  Pa-1s-1).  A, 


VI.  Discussion 

The  goal  of  this  work  was  to  extend  our  source  charac¬ 
terization  work  in  experimental  microwave  inverse  scatter¬ 
ing  [13] — [15]  to  the  problem  of  characterizing  commercial 
transducer  probes  for  transmission  acoustic  inverse  scat¬ 
tering.  The  two  main  challenges  were  1)  the  small  spatial 
scales  of  acoustic  waves  and  2)  the  need  to  characterize 
transducers  using  measurements  instead  of  simulation. 
These  two  differences  motivated  the  derivation  of  a  self¬ 
characterization  method  for  commercial  probes  based  on 
nonlinear  inversion  of  an  acoustic  propagation  model.  We 
derived  the  propagation  model  in  Section  II-C  and  the 
nonlinear  inversion  in  Section  III-C.  The  method  allowed 
us  to  estimate  the  transducer  transmit  coefficients  from 
pair-wise  measurements.  We  further  derived  a  modified 
scattered  field  VIE  in  Section  II-D  so  that,  once  the  trans¬ 
mit  coefficients  were  estimated,  the  incident  fields  com¬ 
puted  with  them  were  rigorously  linked  to  volume  inte¬ 
grals  of  the  inverse  scattering  algorithm.  We  also  used  an 
additional  optimization  with  time  of  flight  to  solve  for  the 
positions  of  the  two  probes.  Several  experimental  assump¬ 
tions  were  necessary  given  the  available  setup,  but  these 
were  limited  as  much  as  possible. 

From  the  results  of  the  characterization  and  imaging, 
together  with  the  assumptions  of  the  experiment,  the 
method  showed  reasonable  success.  The  characterization 
inversion  fit  the  data  well.  The  targets  were  mostly  recov¬ 
ered  and  the  signs  of  their  quantitative  values  followed 
physical  expectations.  Without  the  transducer  model,  the 
inversion  produced  contrasts  with  flipped  signs  and  con¬ 
trasts  that  were  too  large  by  an  order  of  magnitude. 


A  major  assumption  affecting  image  quality  and  inter¬ 
pretation  was  imaging  2-D  slices  of  inherently  3-D  objects. 
This  was  necessary  given  the  limited  data  acquisition  ge¬ 
ometry,  where  we  assumed  the  cylindrical  geometry  and 
long  distances  would  allow  us  to  ignore  out  of  plane  scat¬ 
tering.  This  assumption  likely  contributed  to  the  contrast 
values  being  too  high.  To  illustrate  this,  we  ran  an  ideal¬ 
ized  simulation  showing  the  effects  of  object  reconstruc¬ 
tion  in  a  thin-3-D  domain  when  scattered  fields  are  gener¬ 
ated  from  a  thick-3-D  domain.  Ideal  point-sources  in  the 
same  geometry  as  the  experiments  are  used.  The  object 
consists  of  two  cylinders,  both  with  compressibility  con¬ 
trast  of  0.1,  shown  in  Fig.  14(a).  The  thick-3-D  domain  is 
8  x  8  x  10  A,  whereas  the  thin-3-D  domain  is  8  x  8  x  1  A, 
where  A  =  0.4  mm.  Scattered  fields  are  generated  from  the 
cylinders  in  the  thick  domain.  Fig.  14(b)  is  the  2-D  profile 
of  the  object  reconstructed  in  the  thick-3-D  domain;  Fig. 
14(c)  shows  the  object  reconstructed  in  the  thin-3-D  do¬ 
main.  The  peak  contrast  of  the  object  reconstructed  in  the 
thin  domain  is  overestimated  by  a  factor  of  9.  In  simulated 
examples,  such  as  this,  overestimation  scales  with  the  ra¬ 
tio  of  domain  thicknesses,  because  the  algorithm  is  trying, 
in  effect,  to  conserve  the  scattered  field  power  radiated  by 
the  objects  in  each  domain. 

Finally,  we  did  not  have  independent  confirmation  of 
the  accuracy  of  the  forward-scattering  model,  i.e.,  the 
combination  of  the  estimated  transmit  coefficients,  com¬ 
puted  incident  fields,  and  scattered  field  predictions  by  the 
new  VIE.  This  would  have  required  the  precise  position 
and  material  properties  of  known  targets,  which  was  not 
possible  with  this  experimental  setup.  Motion  undoubt¬ 
edly  contributed  to  image  artifacts  as  well. 


VII.  Conclusion 

We  derived  and  provided  initial  experimental  testing  of 
a  new  self-characterization  method  for  commercial  ultra¬ 
sound  probes  which  was  designed  for  transmission  inverse 
scattering.  We  first  modified  the  traditional  volume  inte¬ 
gral  equations  to  make  them  consistent  with  the  trans¬ 
ducer  model.  This  was  an  extension  of  our  previous  work 
modeling  antennas  for  electromagnetic  inverse  scattering. 
We  then  derived  an  optimization  procedure  based  on  the 
nonlinear  inversion  of  an  acoustic  propagation  model  to 
obtain  the  transducer  transmit  coefficients  from  measured 
data.  We  described  our  two-probe  experiment,  including 
provisions  to  estimate  the  probe  locations  and  align  a  ro¬ 
botic  rotator.  The  characterization  results  showed  that 
the  inversion  fit  the  measured  data  well.  We  then  tested 
our  inverse  scattering  algorithm  with  modified  volume  in¬ 
tegral  equations.  The  images  showed  that  the  values  of 
the  contrasts  of  simple  objects  were  physically  consistent 
with  expected  values,  but  image  quality  suffered  from  as¬ 
sumptions  inherent  in  the  experimental  setup.  This  is  the 
groundwork  for  future  better  versions  of  full- wave  calibra¬ 
tion  with  other  transducers  and  imaging  algorithms,  with- 
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Fig.  14.  Comparison  of  reconstructions  on  domains  of  different  thickness 
given  identical  scattered  fields,  (a)  Object:  two  cylinders  with  compress¬ 
ibility  contrast  of  0.1.  Scattered  fields  are  generated  from  a  domain  with 
a  thickness  of  10A.  (b)  2-D  profile  of  the  reconstruction  on  a  lOA-thick 
domain,  (c)  2-D  profile  of  the  same  object  reconstructed  on  a  lA-thick 
domain;  the  contrasts  are  overestimated  by  a  factor  of  6  (note  an  in¬ 
creased  colorscale).  A 


out  which  we  will  not  be  able  to  improve  ultrasound  much 
beyond  what  has  been  the  state  of  the  art  for  many  years. 

The  transducer  model  and  modified  VIE  are  rigorous 
and  versatile,  but  further  investigation  is  required  to  bet¬ 
ter  evaluate  the  performance  of  the  nonlinear  character¬ 
ization  and  acoustic  imaging  algorithm.  In  future  work, 
we  will  explore  using  simulation  to  obtain  the  transducer 
transmit  coefficients,  and  explore  viable  options  for  in¬ 
dependent  validation  of  the  forward-scattering  model  at 
ultrasound  spatial  scales.  The  propagation  model  could 
also  be  used  to  include  the  full  radiation  and  phase  pat¬ 
terns  of  hydrophones  for  hydrophone  characterization  of 
transducers. 


g3\a,b)  =  aJ-Na^bj,  (61) 


where  a  and  b  are  model  vectors  for  receiver  and  trans¬ 
mitter  positions,  respectively.  After  substituting  the  up¬ 
date  step,  and  rearranging, 


g(mJ 


g(mre_1,m„_1)  -  a„(g(m„_1,v„)  +  g(v„,m„.,)) 
+  <^g(vn,vj. 

(62) 


Substituting  this  into  the  cost  function,  we  have 

2Umn)  =  Ik-1  -  anSl  +  areg 2  Id  +  lCn— 1  “  «nv  Jlk 

(63) 

where 


r„-i  =  g(m„_1,m„_1)  -  d 

(64) 

gl  =  g(mn— L>VJ  +  g(vn>m»-l) 

(65) 

g2  =  g(v„,V„) 

(66) 

C„_!  =  Ill,,  ,  -  111,,. 

(67) 

The  term  g(mn_i,vn)  treats  mn_i  as  receivers  and  vn 
as  sources,  whereas  g(vn,mn_i)  does  the  opposite.  Writing 
out  the  vector  norms 

2T(m„)  =  ||r|£  -  a(r, g^o  +  a2(r,g2)D 

-  a(gikD  +  «2 ||gillb  -  a3(gi,g2)D 

(68) 

+  a2(g2>r)D  -  “3(g2>gl)D  +  ^  4  II  g  2  II D 

+  IIcIIm  -a(c,v)M  -a(v,c)M  +  a2\\vfu 

and  using  the  fact  that  (a,  b)  =  (b,a)*,  we  can  write  this  as 
a  quartic  polynomial  in  an: 

F(x)  =  ax 4  +  bx 3  +  cx 2  +  dx  +  e,  (69) 


Appendix 
Step  Length 

We  choose  an  to  minimize  the  cost  function  at  each 
step.  Substituting  the  expression  for  mn  into  the  cost 
function, 

2i?(m n)  =  llg(m„_i  -  anv n)  -  d|£ 

2  (60) 

"b  ll-^n— 1  n  ^all]V[' 


where 


x  =  a 

(70) 

II 

oT 

to 

5~to 

(71) 

b  =  -2sR(ghg,)D 

(72) 

gill2  +  25?(r,g2)D  +  \\vfM 

(73) 

-  —  25R(r,  gx)D  -  23i(c,v)M 

(74) 

e  =  lr  D  +  IcIm  • 

(75) 

Because  this  forward  model  is  only  polynomial  non¬ 
linear,  we  will  obtain  a  polynomial  in  <a,  the  minimum  of 
which  can  be  found  using  standard  algebraic  techniques. 
To  help  the  bookkeeping,  we  first  write  the  forward  model 
in  matrix  notation  with  two  generic  model  vectors  repre¬ 
senting  receivers  and  transmitters 


Next,  we  must  find  the  minimum  of  F(x)  with  respect 
to  x.  The  first  and  second  derivatives  of  the  quartic  are 
given  by 

F\x)  =  4  ax3  +  3  bx2  +  2  cx  +  d  (76) 

F"(x)  =  12  ax2  +  6bx  +  2c.  (77) 


14 


IEEE  TRANSACTIONS  ON  ULTRASONICS,  FERROELECTRICS,  AND  FREQUENCY  CONTROL,  VOL.  TBC,  NO.  TBC,  TBC  TBC 


The  leading  term  of  the  second  derivative  shows  us  that 
this  quartic  will,  globally,  always  be  convex,  because  a  is 
a  square  vector  norm  and  always  greater  than  zero.  For 
a  convex  quartic,  there  is  either  one  global  minimum,  or 
there  are  at  most  2  local  minima  and  1  local  maximum. 
We  can  separate  these  cases  by  finding  the  zeros  of  F{x), 
for  which  there  is  either  1  real  root  and  two  complex  con¬ 
jugate  roots  (i.e.,  1  global  minimum  of  the  quartic),  or  3 
real  roots  (i.e.,  2  local  minima  and  1  local  maximum  of  the 
quartic).  In  the  latter  case,  we  find  the  global  minimum 
by  back  substitution;  any  ambiguity  that  might  arise  from 
identical  local  minima  (i.e.,  different  possible  step  lengths) 
has  not  been  observed  in  practice  to  be  a  problem.  A  cubic 
is  one  of  the  few  polynomials  for  which  the  zeros  can  be 
determined  analytically.  There  are  several  methods,  which 
can  be  found  in  standard  math  texts. 
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One  application  of  acoustic  droplet  vaporization  (ADV),  a  method  of  converting  biocompatible 
microdroplets  into  microbubbles,  is  to  enhance  locally  high  intensity  focused  ultrasound  (HIFU) 
therapy.  Two  objectives  are  pursued  here:  (1)  the  controlled  creation  of  a  bubble  trench  prior  to 
HIFU  using  ADV  and  (2)  use  of  the  trench  for  increasing  ablation  volumes,  lowering  acoustic 
powers,  and  decreasing  therapy  duration.  Thermally  responsive  phantoms  were  made  with 
perfluorocarbon  emulsion.  Compound  lesions  were  formed  in  a  laboratory  setting  and  a  clinical 
magnetic  resonance  imaging  (MRI) -guided  HIFU  system.  Linear  and  spiral  patterned  compound 
lesions  were  generated  in  trenches.  A  larger  fraction  of  the  HIFU  beam  is  contained  to  increase 
the  generation  of  heat.  Using  the  laboratory  system,  a  90  mm  linear  length  spiral  trench  was 
formed  in  30  s  with  mechanical  beam  steering.  Comparatively,  the  clinical  HIFU  system  formed 
a  19.9mm  linear  length  spiral  trench  in  approximately  Is  with  electronic  beam  steering. 

Lesions  were  imaged  optically  and  with  MRI.  A  uniform  thermal  ablation  volume  of  3.25  mL 
was  achieved  in  55.4  s  (4-times  faster  than  standard  clinical  HIFU  and  14-times  larger  volume 
versus  sum  of  individual  lesions).  Single  lesions  showed  a  400%  volume  increase. 

© 2014  Acoustical  Society  of  America.  [http://dx.doi.org/10.1121/L4828832] 
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I.  INTRODUCTION 

High  intensity  focused  ultrasound  (HIFU)  therapy  is  a 
non-invasive  and  non-ionizing  technique  wherein  focused 
ultrasound  beams  are  emitted  from  a  high-powered  trans¬ 
ducer  to  thermally  ablate  the  tissue  volume  of  interest 
(Crouzet  et  al. ,  2010).  HIFU  systems  commonly  operate  in  a 
frequency  range  of  0.5-5  MHz,  generating  high  focal  inten¬ 
sities  up  to  1-lOkW/cm2  (Shaw  and  Hodnett,  2008).  Such 
intensities  can  increase  tissue  temperature  beyond  60  °C 
within  the  focal  volume  in  seconds,  resulting  in  irreversible 
protein  denaturation,  cell  destruction,  and  associated  coagu- 
lative  necrosis  (Hill  and  ter  Haar,  1995).  Two  principal 
mechanisms,  direct  absorption  of  the  transmitted  pressure 
wave  and  acoustic  cavitation,  can  cause  tissue  heating  syn- 
ergistically  (Coussios  and  Roy,  2008).  The  bioeffects  of 
focused  ultrasound  on  tissues  have  been  investigated  for  dec¬ 
ades  (Bamber  and  Hill,  1979;  Parker,  1983;  Frizzell,  1988) 
and  the  physical  principles  associated  with  HIFU  are  becom¬ 
ing  better  understood  as  HIFU  becomes  an  accepted  modal¬ 
ity  for  treatment  (ter  Haar  and  Coussios,  2007). 
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Focused  HIFU  energy  is  sufficient  to  induce  thermal 
coagulation  in  a  small  volume  within  seconds,  and  therefore 
blood  perfusion  effects  acting  as  thermal  sinks  are  minimal 
(Billard  et  al .,  1990).  Although  promising  and  effective  for 
small  tumors,  attempts  to  expand  the  application  of  HIFU 
thermal  ablation  to  larger  masses  have  had  limited  success 
(Yamakado  et  al .,  2006;  Gervais  et  al .,  2003).  A  major  diffi¬ 
culty  has  been  creating  homogenous,  reproducible,  and  uni¬ 
formly  shaped  lesions  in  a  time  period  that  is  clinically 
manageable  for  treating  large  volumes.  In  addition,  some 
patients  complain  about  local  pain  after  HIFU  therapy, 
which  indicates  normal  tissue  overheating,  especially  at  sites 
of  large  acoustic  impedance  mismatch  (Li  et  al .,  2009; 
Zhang  et  al .,  2010b).  Therefore,  methods  localizing  treat¬ 
ment  and  thus  protecting  sensitive  tissues  would  greatly  ben¬ 
efit  HIFU  and  potentially  yield  new  treatment  locations. 

It  is  well  known  that  the  thermal  contribution  of  micro¬ 
bubbles  can  be  significant.  Ultrasound  contrast  agents 
(UCAs)  have  been  investigated  for  the  enhancement  of  ther¬ 
mal  effects  at  the  HIFU  focus  (e.g.,  Tung  et  al.,  2006). 
However,  UCAs  also  move  the  greatest  heating  position 
from  the  focus  by  up  to  2  cm  for  a  given  transducer  geometry 
and  aberrating  tissues  intervening  the  acoustic  path  as  well 
as  create  surface  lesions  for  concentrations  greater  than 
0.1%  (v/v),  with  some  damage  well  outside  the  lesion  (Tung 
et  al.,  2006).  Alternatively,  acoustic  droplet  vaporization 
(ADV)  is  an  ultrasound-based  method  of  converting  biocom¬ 
patible  microdroplets  (Fabiilli  et  al.,  2013)  into  microbub¬ 
bles.  When  exposed  to  acoustic  pressures  above  the 
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vaporization  threshold,  injected  perfluorocarbon  microdrop¬ 
lets — which  are  similar  in  size  to  UCAs — produce  micro¬ 
bubbles  that  are  more  than  100-fold  larger  in  volume  than 
the  original  microdroplets,  resulting  in  a  local  increase  in  the 
acoustic  absorption.  Hence,  the  increased  heating  with  ADV 
can  be  limited  to  the  acoustic  focus  where  pressures  exceed 
the  ADV  threshold,  allowing  for  lower  source  intensities 
and/or  shorter  exposure  times  for  given  therapeutic  effects  at 
the  focus.  Increased  absorption  by  bubbles  at  the  focus  also 
reduces  risk  to  distal  tissues  by  decreasing  the  transmitted 
acoustic  energy  (Lo  et  al. ,  2006). 

ADV  bubbles  are  stable  for  several  minutes  or  longer 
and  can  cause  vascular  occlusion  when  lodged  in  a  capillary 
or  other  small  vessels  (Kripfgans  et  al .,  2000),  which  is 
favorable  for  HIFU  therapy.  These  microbubbles  are  highly 
echogenic  (Zhang  et  al .,  2010a,  Kripfgans,  2002),  thus  facili¬ 
tating  real-time  image-guided  therapy.  Therefore,  the  ultra¬ 
sonic  monitoring  of  microbubbles  created  during  the  ADV 
process  may  offer  opportunities  to  visualize  the  location  and 
spatial  extend  of  the  HIFU  process.  Additionally,  a  dense 
bubble  plane  forms  a  barrier  for  acoustic  waves  and  impedes 
wave  propagation  to  distal  tissues  and  can  double  pressure 
amplitudes  proximal  to  the  bubble  wall  (Lo  et  al .,  2006). 
Similar  studies  by  Zhang  et  al.  (2013)  demonstrated  in  vitro 
that  the  inclusion  of  0.008%  (v/v)  phase-shift  nanoemulsion 
(PSNE)  in  phantoms  significantly  reduced  the  power 
required  to  generate  9  mm3  lesions  compared  to  HIFU  with¬ 
out  PSNE. 

This  study  is  focused  on  substantially  enhancing  HIFU 
by  increasing  thermal  lesion  volume  (or  equivalently  short¬ 
ening  treatment  time)  through  precisely  controlled  produc¬ 
tion  of  ADV  microbubbles  from  micrometer- sized  droplets. 
This  technology  is  expected  to  be  applicable  across  a  broad 
range  of  current  and  future  uses  of  HIFU  ablative  therapy. 

II.  MATERIALS  AND  METHODS 
A.  Thermal  phantom  preparation 

Thermally  responsive  phantoms  were  fabricated  based 
on  previously  described  methods  (Takegami  et  al .,  2004). 
The  phantoms  consisted  of  33%  (v/v)  of  an  aqueous  acryl¬ 
amide  solution  (30%  wt/vol,  Sigma-Aldrich,  St.  Louis,  MO), 
31.4%  (v/v)  of  degassed  deionized  water,  35%  (v/v)  of  egg 
white,  0.5%  (v/v)  of  10%  ammonium  persulfate  (Sigma- 
Aldrich),  and  0.1%  (v/v)  tetramethylethylenediamine 

(Sigma-Aldrich).  Lipid  droplets  were  made  according  to  a 
previously  described  protocol  (Fabiilli  et  al .,  2009)  and 
stored  in  the  refrigerator  (5  °C)  for  up  to  4  weeks.  Prior  to 
use,  the  droplet  size  distribution  and  number  density  were 
measured  with  a  Coulter  counter  (Multisizer  3,  Beckman 
Coulter  Inc.,  Fullerton,  CA).  Approximately  99%  of  droplets 
in  the  emulsion  were  smaller  than  8  pm  in  diameter.  The 
mean  diameter  of  the  droplets  was  2.0  ±0.1  /un. 

During  phantom  preparation,  the  gel  solution  was 
degassed  before  the  addition  of  droplets  to  minimize  spurious 
gas  inclusions  that  could  act  as  cavitation  nuclei.  The  concen¬ 
tration  of  droplets  within  the  phantom  was  105  droplets/mL, 
which  corresponds  to  0.00004%  (v/v).  These  egg  white-based 
phantoms  were  optically  transparent.  If  HIFU  increased  the 


temperature  in  a  portion  of  the  phantom  to  >60  °C  for  >5  s, 
the  egg  white  protein  would  denature  and  coagulate,  resulting 
in  a  permanently  opaque,  optically  observable  lesion.  The 
acoustic  attenuation  of  these  phantoms  was  approximately 
0.3  dB/cm  at  the  resonance  frequency  of  our  HIFU  transducer 
(1.44  MHz),  which  is  also  very  close  to  the  frequency  of  the 
clinical  HIFU  system  of  1 .45  MHz. 


B.  Laboratory  HIFU  lesion  generation 

The  HIFU  therapy  system  consisted  of  a  high  power 
spherical  section  transducer  (63.5  mm  diameter,  f-number  of 
1,  resonance  frequency  1.44  MHz,  Etalon  940501,  Lebanon, 
IN)  driven  with  a  power  amplifier  (A-300,  ENI  Inc., 
Rochester,  NY).  Low  power  pulse-echo  operation  of  the 
HIFU  transducer  allowed  for  focal  localization  in  the  phan¬ 
tom  by  detection  of  phantom  boundaries  by  use  of  a  micro 
positioning  system.  Also  under  pulse-echo  guidance,  a  nee¬ 
dle  thermocouple  (300  /mi  diameter,  Omega,  Stamford,  CT) 
was  inserted  through  a  guide  into  the  side  of  the  phantom 
2  mm  away  from  the  transducer  focus,  and  a  data  logger 
(TC-08,  Pico  Technology,  Cambridge,  UK)  was  used  to  re¬ 
cord  the  temperature.  A  function  generator  (3314A,  Agilent, 
Palo  Alto,  CA)  was  gated  with  a  secondary  function  genera¬ 
tor  (33120A,  Agilent)  to  produce  continuous  wave  (CW)  sig¬ 
nals  to  the  amplifier  for  ultrasound  exposures  with  durations 
of  5  s.  The  lateral  beam  width  and  axial  depth  of  field  at  the 
focus  [full-width  at  half  maximum  (FWHM)  pressure]  were 
1.6  mm  and  8.2  mm,  respectively.  The  therapy  transducer 
was  calibrated  using  a  membrane  hydrophone  (ST223,  Sonic 
Industries,  Hatboro,  PA)  in  degassed  water  at  room  tempera¬ 
ture,  and  the  output  waveforms  from  the  hydrophone  were 
digitized  (9314L,  LeCroy  Corp.,  Chestnut  Ridge,  NY).  A 
peak  negative  pressure  of  7.4  MPa  and  a  peak  positive 
pressure  of  17.7  MPa  were  measured  at  the  focus.  The 
maximum  focal  intensity  was  then  estimated  from  the 
measured  peak  positive  pressure  (Shaw  and  Hodnett, 
2008)  by  approximating  the  2D  acoustic  beam  profile 
with  a  Gaussian  function  and  assuming  linear  acoustic 
(first  order  assumptions) 


p(x,y)  =  p+  e-(x  /2<7 *  +y  /2^  ), 


p(x,y)2dx  dy 


/Sptp 


n(FWHM  /  2) 


2  ’ 


where  p(x,y)  is  the  focal  plane  pressure  distribution,  p+  is  the 
peak  positive  pressure  amplitude  of  17.7  MPa,  ox  and  oy  are 
related  to  the  FWHM  beam  width  by  =  FWHM/2.35, 
^acoustic  is  the  acoustic  power  to  be  determined,  Z  is  the 
wave  impedance  of  1.5  MRayl,  dx  times  dy  compose  an  area 
element  dA  in  the  focal  plane  and  /sptp  is  the  spatial  peak, 
temporal  peak,  acoustic  intensity  derived  from  the  total 
acoustic  power  over  the  FWHM  focal  plane  area.  For  the 
laboratory  HIFU  system,  Pac oustic  is  computed  to  304  W  and 
/sptp  to  15.1  kW/cm  . 
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(A)  Backwall  creation  (B)  Trench  generation 


16  mm 

(C)  Individual  HIFU  exposure  (D)  Resulting  composite  lesion 


FIG.  1 .  (Color  online)  Three-dimensional  schematic  depiction  of  the  acous¬ 
tic  trench  created  for  HIFU.  (A)  First  a  back  wall  bubble  layer  is  created  by 
means  of  Acoustic  Droplet  Vaporization.  (B)  Then  a  bubble  spiral  wall  is 
created  on  top  of  it  to  form  a  trench.  (C)  Inside  the  bubble  trench,  the  acous¬ 
tic  field  of  a  HIFU  beam  is  used  to  create  thermal  lesions.  Discrete  lesions 
are  placed  along  the  spiral  track  of  the  trench  and  result  in  a  solid  lesion  of 
the  shape  of  a  thick  disk  of  15  mm  height  (D).  Trench  geometries  depend  on 
the  HIFU  point  spread  function  and  steering  capabilities  of  the  HIFU  beam. 
Note:  Approximate  aspect  ratios  are  shown.  The  beam  widths  of  the  clinical 
and  laboratory  HIFU  systems  are  1.0  and  1.6  mm,  respectively,  and  the 
trench  walls  are  spaced  approximately  5  mm.  A  3.25  mL  lesion  can  be  cre¬ 
ated  in  55.4  s.  Actual  heating  required  40  s  while  the  creation  of  the  acoustic 
trench  required  15.4  s. 

A  programmable  3-axis  positioning  system  was  used  to 
translate  the  focused  transducer  in  a  spiral  pattern  as  illus¬ 
trated  in  Fig.  1.  A  short  tone  burst  of  10  cycles  at  1.44  MHz 
with  a  0.25%  duty  cycle  was  applied  to  the  therapy  trans¬ 
ducer  to  vaporize  the  droplets  in  thermal  phantom  blocks 
(5.5  cm  x  5.5  cm  x  6.5  cm),  forming  a  bubble  back  wall 
[10  mm  in  thickness,  Fig.  1(A)]  and  then  a  spiral  sidewall  [3 
mm  in  thickness,  Fig.  1(B)].  The  purpose  of  the  back  and 
sidewalls  is  to  reflect  and  largely  contain  the  fraction  of  the 
incoming  and  bubble- scattered  HIFU  beam,  confining  a 
greatly  increased  fraction  of  the  heat  generation  within  the 
trench. 

The  formation  of  the  spiral  shaped  trench  was  achieved 
within  30  s.  Subsequently,  compound  lesions  were  formed 
by  5  s  exposures  of  the  HIFU  beam  inside  the  trench  [Fig. 
1(C)]  separated  by  5  s  delays.  As  shown  here,  ultrasound 
was  focused  at  discrete  locations  along  the  trench  thereby 
generating  a  solid  thermal  lesion  disk  [Fig.  1(D)].  Center-to- 
center  spacing  between  neighboring  lesions  within  the  trench 
was  5.5  mm  (right  column  in  Fig.  2).  The  geometry  of  the 
acoustic  trenches  was  verified  prior  to  HIFU  exposure,  using 
a  GE  Logiq  9  ultrasound  scanner  (GE  Ultrasound, 
Milwaukee,  WI)  with  a  7  L  array  operating  at  7  MHz. 

C.  Lesion  size  estimation  using  visible  macroscopic 
imaging  and  MRI 

All  lesions  were  imaged  at  2  T  with  T2- weighted  magnetic 
resonance  imaging  (MRI)  with  1  mm  x  0.5  mm  x  0.25  mm 
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resolution,  and  optical  macroscopic  observations  in  sliced 
phantoms.  In  MRI,  the  transverse  magnetization  relaxation 
time  constant,  T2,  is  sensitive  to  slow  molecular  motions  of 
water.  The  denatured  and  cross-linked  proteins  caused  a  sig¬ 
nificant  decrease  in  the  water  proton  mobility  and  thereby  a 
decrease  in  the  water  proton  T2  (Lambelet  et  al. ,  1991), 
which  was  observed  as  a  darker  region  in  Fig.  2.  Therefore, 
MRI  provided  an  alternative,  independent,  and  direct  method 
of  measuring  protein  coagulation,  and  it  was  utilized  to 
assess  the  uniformity  of  compound  thermal  lesions.  The 
spin-echo  MRI  imaging  sequence  with  T2  weighting  was 
performed  within  1  h  after  thermal  lesion  production  in  the 
phantoms. 

The  repetition  time  TR,  the  amount  of  time  between 
successive  pulse  sequences  applied  to  the  same  slice,  was 

4.5  s.  The  echo  time  TE,  the  time  between  the  90°  pulse  and 
the  peak  of  the  spin  echo  signal  and  inversion  recovery 
pulse,  was  120  ms.  After  MRI,  the  lesions  were  sliced  and 
dissected  based  on  the  visible  color  change.  Lesion  volumes 
of  the  laboratory  HIFU  system  were  also  estimated  by  a  fluid 
displacement  method. 

D.  Lesion  formation  using  an  MRI  guided  clinical  HIFU 
system 

A  Philips  Sonalleve  MRI-HIFU  system  utilizing  a  256- 
element  phased  array  was  used  to  generate  HIFU  lesions. 
For  the  generation  of  spiral  shaped  trenches  a  first  generation 
Sonalleve  VI  system  with  a  120  mm  HIFU  transducer  focal 
length,  implemented  on  a  1.5  T  MR  scanner,  was  employed. 
A  newer  system,  a  Sonalleve  V2  (3T  MR  scanner,  140  mm 
HIFU  transducer  focal  length),  was  employed  for  the  genera¬ 
tion  of  single  lesions.  Both  Sonalleve  systems,  while  not 
identical,  are  capable  of  generating  either  lesion  shape.  At 
the  employed  transmit  center  frequency  of  1 .45  MHz,  a  max¬ 
imum  electronic  steering  of  ±10  mm  was  achieved  at  a 
depth  of  10  cm  for  spiral  pattern  generation  analogous  to  that 
of  the  single  element  transducer.  Emulsions  were  vaporized 
by  use  of  continuous  wave  excitation.  For  the  creation  of  the 
bubble  back  wall  and  bubble  spiral  trench  using  the 
Sonalleve  VI  system,  an  exposure  of  50  W  was  employed 
for  100  ms.  For  the  creation  and  characterization  of  single 
thermal  lesions,  the  Sonalleve  V2  system  was  employed 
using  exposures  ranging  from  30  to  200  W  for  20  s  at 
1.45  MHz  center  frequency.  Calibrations  under  free  field 
conditions  and  1  W  excitation  in  a  water  tank  yielded  a  pres¬ 
sure  amplitude  of  1.12  MPa  and  a  point  spread  function 
FWHM  intensity  of  1.00  mm  lateral  and  7.99  mm  axial. 
Under  the  assumption  of  linear  propagation,  200  W  excita¬ 
tion  in  water  would  result  in  a  pressure  amplitude  of 
15.9  MPa.  Assuming  that  66%  of  the  beam  intersects  the 
FWHM  of  the  focus  as  a  planar  wave  traveling  in  water  with 

1.5  MRayl  impedance,  this  pressure  corresponds  to 
1 2.7  kW/cm2(  16.2  MPa). 

Lesion  volumes  were  measured  by  MRI  and  computed 
as  all  volumetric  pixels  with  T2  values  deviating  from  back¬ 
ground  as  described  by  Zhang  et  al.  (2011).  Lesion  volume 
increases  due  to  the  presence  of  perfluorocarbon  emulsions 
in  the  phantoms  during  HIFU  were  quantified  as  the  ratios  of 
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HIFU  with  ADV 


HiFU  without  ADV 


FIG.  2.  (Color  online)  Optical  (top 
row)  and  T2-weighted  MR  images 
(bottom  two  rows)  of  lesions  (i.e., 
denatured  protein)  generated  within 
thermal  phantoms  using  the  laboratory 
HIFU  system.  The  left  and  right  col¬ 
umns  depict  lesions  generated  with  and 
without  ADV-enhancement  (see  Fig. 
1),  respectively.  The  scale  bar  is  1  cm. 
The  solid,  double  and  dashed  lines 
superimposed  to  the  middle  and  bot¬ 
tom  row  MRI  scans  indicate  where  the 
lateral-elevational  and  the  lateral-axial 
images  intersect. 


lesion  volume  with  and  without  emulsions  present.  Lesion 
sizes  were  quantified  as  a  function  of  applied  acoustic  power 
ranging  from  30  to  200  Watts  and  in  situ  temperature  moni¬ 
toring  with  a  therapeutic  endpoint  of  75  °C  (see  Fig.  3). 
Statistically  significant  differences  between  lesions  sizes 
(i.e.,  with  and  without  ADV)  were  determined  using  a 
Student’s  Ftest.  A  significance  level  of  0.05  was  used  for  all 
comparisons. 

In  Fig.  3(A),  real-time  heating  is  assessed  and  then  dis¬ 
played  via  a  pseudo  color  map  whose  scale  is  shown  on  the 
right  colorbar.  The  temperature  map  is  overlaid  on  the  con¬ 
tinuously  updated  magnitude  image  resulting  from  the  gradi¬ 
ent  echo  sequence  that  is  used  to  measure  temperature.  Since 


the  thermal  map  scan  sequence  is  not  optimized  for  good 
contrast,  but  rather  for  fast  acquisition  and  good  signal  to 
noise  at  the  center,  the  image  appears  gray  and  is  subject  to 
some  zebra  artifacts  (Stadler  et  al .,  2007)  at  the  edge  of 
image.  In  this  sequence,  ADV  was  performed  at  an  output  of 
100  W.  A  temperature  of  83  °C  was  reached  in  part  of  the 
lesion.  Conversely,  Fig.  3(B)  displays  the  real-time  tempera¬ 
ture  mapping  employing  a  static  T2  weighted  spin  echo 
image  sequence  acquired  before  the  heating  sequence.  The 
T2  weighted  spin  sequence  is  slow  but  provides  better  con¬ 
trast  without  artifact  or  image  distortion,  which  allows  dis¬ 
crimination  of  the  various  layers  in  the  phantom  and  its 
surroundings  more  clearly. 


(A)  {8} 

FIG.  3.  (Color  online)  MRI-guided  temperature  readout  during  HIFU  exposure.  (A)  A  temperature  map,  scale  bar  on  the  right  side,  is  overlaid  on  the  magni¬ 
tude  image  of  the  MRI  gradient  echo  sequence.  The  gradient  echo  sequence  is  optimized  for  fast  acquisition  and  signal-to-noise  at  the  center.  (B)  Conversely, 
using  a  T2  weighted  spin  sequence  once,  a  static  background  image  is  generated  and  placed  behind  the  continuously  updated  heating  color  map,  therefore  pro¬ 
viding  better  contrast  for  distinguishing  phantom  details  and  other  important  structures.  Note:  The  HIFU  beam  enters  from  the  left  into  the  phantom  and  the 
spatial  extent  of  the  beam  is  illustrated  by  the  two  shaded  opposing  triangles. 
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HIFU  with  ADV 


HIFU  without  ADV 


FIG.  4.  (Color  online)  Optical  images  (side  view)  of  linear  compound 
lesions  (white)  in  doped  phantoms,  generated  with  the  laboratory  system, 
with  (left)  and  without  (right)  a  bubble  wall.  The  bubble  back  wall  is  visible 
below  the  lesion  on  the  left.  Note  that  the  tadpole  shape  characteristic  of 
bubble-enhanced  HIFU  lesions  is  not  evident  with  the  bubble  wall.  A  ther¬ 
mal  ablation  volume  of  15  mL  was  achieved  in  4  min  with  ADV  and  the 
acoustic  trench.  The  scale  bar  is  1  cm.  The  dashed  arrowhead  line  shows 
propagation  direction  of  the  HIFU. 

III.  RESULTS  AND  DISCUSSION 
A.  Laboratory  HIFU  lesion  generation 

Compound  lesions  in  a  spiral  pattern  were  formed  in 
phantoms  with  and  without  droplets,  in  the  following 
referred  to  as  “doped  phantom”  and  “control  phantom,” 
respectively.  Initial  results,  shown  in  Fig.  2,  were  similar  for 
both  optical  macroscopic  and  T2- weighted  MR  images.  In 
MR  images,  a  high  contrast  between  the  mean  gray  levels 
inside  the  lesions  and  the  surrounding  background  was 
observed,  which  showed  a  difference  of  at  least  nine  stand¬ 
ard  deviations.  The  coefficients  of  variance  were  found  to  be 
similar  in  the  lesion  and  in  the  background  medium  (i.e.,  the 
standard  deviations  did  not  exceed  5%  of  the  mean  value). 
Therefore,  a  simple  amplitude  thresholding  was  applied  in 
each  section  to  discriminate  pixels  included  in  the  lesions 
from  pixels  of  the  background,  with  a  threshold  set  to  the 
mid-value  between  the  mean  gray  level  values  in  the  lesion 
and  in  the  background. 

For  lesions  generated  in  a  spiral  pattern  (Fig.  2),  a  2.5- 
fold  increase  in  the  axial  dimension  of  the  lesion  was 
obtained  with  the  acoustic  trench.  Figure  4  shows  optical 
images  of  linear  compound  lesions.  Given  the  homogene¬ 
ously  complete  filling  of  the  treatment  space,  there  is  a  con¬ 
siderably  larger  volume  that  is  treated  when  using  the  ADV 
enhancement  and  acoustic  trench  approach.  For  instance, 
here  a  uniform  thermal  ablation  volume  of  15  mL  was 


achieved  in  4  min  with  ADV  and  the  acoustic  trench.  A  com¬ 
mon  term  in  the  HIFU  literature  is  the  time  spent  on  a  lesion 
of  a  given  volume.  Here  4  min  for  a  15  mL  lesion  yields 
16  s/mL. 

In  both  optical  and  MR  images  of  the  HIFU  lesion  in 
the  presence  of  ADV  bubbles,  the  back  bubble  wall  was  visi¬ 
ble  below  the  lesion.  The  tadpole  shape  characteristic  of 
bubble-enhanced  HIFU  lesions  was  not  evident  with  the  bub¬ 
ble  wall  in  place  (Fig.  4).  As  seen  in  Fig.  4,  the  thermal 
lesion,  indicated  by  the  opaque  denatured  egg-white,  ended 
abruptly  and  uniformly  across  the  back  wall.  The  layer  of 
this  exposure  plane  (i.e.,  proximal  to  the  transducer)  can 
then  be  the  back  wall  for  the  next  layer.  It  is  expected  that 
ADV  occurs  during  lesion  formation  as  both  ADV  and  lesion 
formation  were  executed  with  the  same  pressure  amplitudes 
though  longer  exposure  times  were  required  for  lesion  for¬ 
mation.  The  function  of  the  acoustic  trench  sidewalls  is  to 
reflect  the  lateral  and  elevational  beam  components  hereby 
resulting  in  a  homogeneous  lesion  versus  the  tadpole  shape. 

B.  Lesion  formation  results  using  a  clinical  HIFU 
system 

Rapid  repetition  of  electronically  steered  therapy  pulses 
(40  pulses  per  second,  50  W)  in  the  clinical  HIFU  system 
allowed  for  the  generation  of  nearly  homogeneous  composite 
lesions  at  17.0  s/mL  or  better  (Fig.  5).  The  volume  of  com¬ 
pound  lesion  (3.24  mL)  is  14-times  larger  than  the  sum  total 
single  volumes  (8  x  28.4  /iL  =  227  /iL).  The  time  to  create  a 
spiral  compound  lesion  is  55  s  (i.e.,  17.0  s/mL),  versus  40  s 
(i.e.,  176s/mL)  to  create  eight  individual  lesions  with  no 
back  wall  and  no  spiral  trench.  Thus,  a  1400%  increase  in 
lesion  volume  was  achieved  with  only  a  37%  increase  in 
treatment  time. 

Consistent  with  earlier  results  on  the  laboratory  system 
(see  also  Zhang  et  al. ,  2011),  for  acoustic  power  levels  rang¬ 
ing  from  30  to  200  W,  lesion  volumes  of  doped  phantoms 
increased  compared  to  lesion  volumes  in  control  phantoms 
(see  Fig.  6).  Lesion  temperature  was  recorded  in  real-time 
and  used  as  the  therapy  end-point  while  insonifying  at  the 
employed  acoustic  power.  Lesion  size  was  measured  over  a 
range  of  power  levels  and  axial  as  well  as  lateral  lesion 
dimensions.  The  minimum  power  level  at  which  lesions 


FIG.  5.  B-mode  images  of  single  and  compound  lesions  from  the  clinical  HIFU  system  are  shown  in  panels  (A),  (B),  and  (C),  respectively.  (A)  Eight  individ¬ 
ual  bubble  clouds  and  HIFU  from  ADV  produced  with  a  5  s  per  cloud  sweep  speed  and  illustrate  the  extent  of  single  lesions.  (B)  A  spiral  shaped  trench  wall  is 
created  within  4  s  by  rapid  succession  of  individual  acoustic  electronically  steered  pulses  for  ADV.  Forty-four  individual  beam  positions  were  targeted  with 
100  ms  each  to  form  the  shown  spiral.  (C)  In  a  second  step,  as  schematically  illustrated  in  Fig.  1,  the  formed  acoustic  trench  was  filled-in  by  eight  additional 
individually  positioned  HIFU  beams  (50  W,  5  s  each),  which  then  created  the  solid  lesion  that  is  shown.  A  backside  wall  (not  shown)  was  created  before  the 
spiral  with  1 10  beams  in  1 1  s.  The  total  lesion  (3.25  mL)  creation  time  including  bubble  clouds  was  55  s  (see  Fig.  1)  or  17.0  s/mL. 
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FIG.  6.  (Color  online)  (A)  Lesion  sizes  were  estimated  using  T2-weighted 
MR  imaging.  Doped  phantoms  (open  circles)  showed  larger  lesion  volumes 
than  control  phantoms  (stars).  To  reach  an  average  lesion  volume  of 
0.105  mL  in  the  control  phantom,  the  acoustic  power  must  be  more  than 
doubled  (213%).  At  150  W  acoustic  power  the  lesion  volume  ratio  of  doped 
versus  control  exceeds  8-times.  Vertical  lines  indicate  the  thresholds,  below 
which  lesions  are  not  observed.  (B)  and  (C)  Axial  and  lateral  lesion  diame¬ 
ters  were  measured.  Lesion  volumes  in  doped  phantoms  showed  a  greater 
functional  gradient  with  applied  acoustic  power  compared  to  control  phan¬ 
toms.  Data  are  shown  as  the  mean  ±  standard  deviation  for  n  =  4.  “#”  indi¬ 
cate  p  <  0.05  versus  control  phantoms. 

were  observed  was  40  W  (doped  phantom),  which  corre¬ 
sponded  to  a  lesion  volume  of  17  ±  11  /iL.  Figure  6  shows 
that  at  150  W  acoustic  power,  doped  phantoms  yielded 
eight-times  larger  lesions  for  the  same  acoustic  power  than 
control  phantoms.  To  reach  an  average  lesion  volume  of 
0.105  mL  in  the  control  phantom,  the  acoustic  power  must 
be  more  than  doubled  (213%)  compared  to  the  power  used 
in  the  doped  phantom.  Doped  lesions  with  volumes  statisti¬ 
cally  different  from  controls  are  marked  with  a  “#”  symbol 
in  Fig.  6(A).  Regression  analyses  [Figs.  6(B)  and  6(C)]  show 
that  both  the  lateral  and  axial  lesion  lengths  increase  with 
^56.0%  greater  slope  when  comparing  doped  to  control 
phantoms.  A  constant  offset  of  the  fitting  function  over  the 
tested  acoustic  power  possibly  indicates  the  larger  depth  of 
field  compared  to  the  lateral  beam  diameter. 

As  shown  in  Fig.  7,  the  lesion  shapes  within  control  and 
doped  phantoms  differed,  with  lesions  within  control  phan¬ 
toms  showing  higher  symmetry  along  the  axial  direction. 
Doped  phantoms  exhibited  lesions  that  progressed  toward 
the  sound  source  and  were  cone  shaped,  except  for  a  tail  dis¬ 
tal  to  the  triangle.  At  each  axial  position  of  a  lesion  the 


lateral-elevational  cross-sectional  area  A  was  measured. 
Then  an  effective  radius  was  calculated  as  reff  =  yjAjn. 
Line  plots  of  reff  are  shown  as  a  function  of  axial  position  to 
illustrate  the  three-dimensional  character  of  the  lesions. 

The  acceleration  of  HIFU  procedures  by  use  of  adminis¬ 
trated  or  induced  gas  bodies  has  been  investigated  by  numer¬ 
ous  research  groups.  Peng  et  al.  (2012)  compared 
ultrasound-guided  HIFU  treatment  for  uterine  fibroids  to 
treatment  enhanced  by  SonoVue,  administrated  10  min 
before  HIFU  as  a  bolus  of  2.0  mL  normal  saline  with 
10.0  mg  SonoVue  via  the  hand  vein.  For  fibroids  smaller 
than  4  cm  in  diameter,  average  sonication  times  decreased 
from  69  s/mL  (HIFU  only)  to  34s/mL  (HIFU  +  SonoVue), 
meanwhile  the  range  decreased  as  well  2-619  s/mL  versus 
15-334  s/mL,  but  is  not  significantly  different  at  p  =  0.354. 
For  ADV  accelerated  HIFU,  in  comparison,  single  lesion 
formation  was  achieved  in  the  clinical  HIFU  system  in  less 
than  20  s,  resulting  in  a  0.67  mL  volume,  hence  allowing  for 
approximately  30  s/mL  lesioning  rates.  For  lesions  based  on 
spiral  trenches,  the  rate  improved  to  17.0  s/mL. 

During  image-guided  ultrasound  therapy,  the  rate  of 
massive  gray  scale  changes  during  lesion  formation  was 
50.0%  (HIFU  only)  versus  77.8%  (HIFU  +  SonoVue)  (Peng 
et  al.,  2012).  Chung  et  al.  (2012)  reported  that  the  volume  of 
coagulative  necrosis  was  significantly  larger  for  contrast- 
assisted  HIFU  in  comparison  to  HIFU  alone.  Here  the  vol¬ 
ume  increased  from  0.83  mL  to  1.75  mL  (p<0.05). 
Simultaneously,  tissue  ablation  times  significantly  decreased 
from  12.1  s  to  8.6  s  (p<  0.05).  Illing  et  al.  (2005)  stated  that 
observed  adverse  events  were  localized  to  the  area  of  treat¬ 
ment.  Following  HIFU  treatment,  skin  toxicity  occurred  in 
27%  of  the  cases  (8/30).  These  were  not  clinically  relevant 
(except  for  one  case)  and  resolved  spontaneously.  Ablation 
at  large  depth  (e.g.,  kidney  tumors)  may  be  more  difficult 
and  less  reliable  than  at  more  superficial  depths  (e.g.,  liver 
tumors).  Illing  et  al.  speculated  that  this  decreased  reliability 
“may  be  due  to  the  greater  depth  of  renal  tumo[u]rs  and  the 
presence  of  a  perinephric  fat  layer,  both  of  which  attenuate 
the  ultrasonic  beam.” 

To  treat  rapidly  enough  for  efficient  use  of  expensive 
facilities  and  staff  as  well  as  reasonable  duration  of  anesthesia, 
it  is  essential  to  enhance  heating  without  damaging  the  skin 
and  other  overlying  tissues.  For  example,  a  single  uterine  fi¬ 
broid  with  a  volume  of  17  mL  requires  approximately  20  min 
of  treatment  time  using  current  standard  HIFU.  Using  ADV 
enhanced  trenches,  this  time  would  decrease  to  5  min.  The 
treatment  time  for  a  6  cm  diameter  hepatocellular  adenoma, 
which  currently  requires  2  h,  would  be  reduced  to  30  min. 

ADV  shows  the  potential  to  provide  local  control  of  heat 
deposition  while  minimizing  tissue  damage  caused  by  micro¬ 
bubbles  in  the  path  of  therapeutic  ultrasound  beams,  as  is  the 
case  in  use  of  UCAs.  In  this  study,  an  appropriate  bubble  den¬ 
sity  was  achieved  consistent  with  therapeutic  ultrasound  while 
providing  uniform  ablation.  The  required  droplet  concentra¬ 
tion  was  comparable  to  diagnostic  UCA  use.  While  the  emul¬ 
sion  core  material,  dodecafluoropentane  (C5F12),  is  an  inert 
fluid  with  high  biocompatibility,  in  vivo  administration  of 
large  quantities  may  cause  bioeffects  (Zhang  et  al.,  2011).  In 
this  paper,  105  droplets/mL  were  used  to  enhance  HIFU 
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FIG.  7.  (Color  online)  (Left)  MRI 
images  showing  coronal  cross-sections 
of  generated  HIFU  lesions  at  200  W 
acoustic  power  and  20  s  exposure  for 
doped  and  control  phantoms.  The 
propagation  direction  of  the  HIFU  was 
from  left  to  right.  Achieved  lesion  vol¬ 
umes  and  volumetric  creation  speeds 
were  0.654  mL  and  30.6  s/mL  as  well 
as  0.185  mL  and  108s/mL  for  doped 
(top)  and  control  (bottom)  phantoms, 
respectively.  (Right)  For  each  lesion 
the  lateral-elevational  cross-section 
was  fit  to  a  circle.  The  radius  of  the  fit¬ 
ted  circle  is  plotted  as  a  function  of 
axial  position  through  the  lesion. 
Lesion  geometry  asymmetries  between 
doped  and  control  phantoms  can  be 
seen. 


lesions.  In  a  highly  perfused  organ  such  as  liver,  with  an  aver¬ 
age  fractional  blood  volume  of  31%  (Schwickert  et  al. ,  2005), 
3.23  x  105  droplets/mL  of  blood  are  required  for  HIFU  lesions 
in  a  70  kg  human  with  5L  blood  volume.  Systemically,  this 
yields  2.3  x  107  droplets/kg,  which  is  282-times  lower  in  con¬ 
centration  than  a  recent  biodistribution  study  where  a  similar 
emulsion  was  safely  tolerated  at  6.5  x  109  droplets/kg 
(Fabiilli  et  al .,  2013). 

Our  phantom  experiments  demonstrated  no  heating 
effects  such  as  protein  denaturation  in  the  pre-focal  or 
post-focal  ultrasound  propagation  pathway  (see  Fig.  3).  In 
addition,  use  of  an  ADV-bubble  back  wall  limits  distal  ultra¬ 
sound  propagation,  thus  limiting  damage  beyond  the  wall.  It 
is  possible  to  create  these  distal  barriers  prior  to  treatment  as 
an  additional  layer  of  protection  for  sensitive  tissues  such  as 
nerves  in  clinical  implementations.  It  is  therefore  hypothe¬ 
sized  that,  building  upon  previously  produced  lesions  with 
ADV  bubbles  should  allow  for  the  reduction  of  the  transmit¬ 
ted  acoustic  energy,  which  can  be  advantageous  in  sensitive 
tissue  regions.  These  bubble  walls  can  be  created  using  low 
duty  cycle,  microsecond  pulses  that  will  not  generate  ther¬ 
mal  lesions.  Further  investigations  of  the  requisite  bubble 
wall  characteristics  and  the  level  of  protection  afforded  by 
the  wall  are  warranted  as  are  studies  of  the  optimal  location 
of  the  ultrasound  beam  focus  for  heating  in  front  of  the  bar¬ 
rier  wall.  Our  results  also  indicated  that  for  the  treatment  of 
a  large  volume,  the  ADV-enhanced  HIFU  would  require 
fewer  targeting  sites  to  create  a  compound  lesion  throughout 
a  predetermined  volume,  because  the  size  of  the  individual 
lesions  was  increased.  Lesion  size  diminishes  with  decreas¬ 
ing  acoustic  power,  as  seen  when  treating  tissues  at  large 
depth.  Compensatory  increase  in  acoustic  power  is  limited 
due  to  its  consequential  damage  of  overlying  tissues. 
However,  an  increase  in  the  maximum  depth  of  treatment  is 
achievable  by  means  of  ADV  enhanced  HIFU  as  the 
required  acoustic  power  is  on  average  lowered  by  47%  for 
the  same  depth  lesion  of  0. 1  mL  volume,  when  compared  to 
non- ADV  HIFU  (Fig.  6).  Therefore,  currently  unbeatable 
regions  due  to  depth  or  restricted  aperture  should  be  realiz¬ 
able  by  the  use  of  ADV. 
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Another  parameter  warranting  further  study  is  the  op¬ 
timum  size  for  the  emulsion.  In  this  paper,  a  micron-sized 
emulsion,  similar  in  size  to  UCA,  was  used.  A  previous 
in  vitro  study  (Zhang  et  al .,  2013)  utilized  a  50-fold 
higher  concentration  of  submicron- sized  emulsion,  which 
could  ultimately  accumulate  within  a  tumor  via  extravasa¬ 
tion  through  leaky  vasculature.  This  type  of  passive  tar¬ 
geting  could  further  enhance  localized  heating  effects, 
though  the  relationship  between  emulsion  size  and  the 
concentration  required  for  enhanced  heating  is  not  fully 
elucidated. 

IV.  SUMMARY  AND  CONCLUSION 

ADV  greatly  enhanced  HIFU  therapy  by  increasing  the 
efficiency  of  focal  energy  deposition.  Large,  compound 
lesions  can  be  generated  at  a  rate  of  17  s/mL  of  lesion  volume. 
Uniform  stacking  of  spiral  lesions  occurred  without  the  distal 
gaps  characteristic  of  lesions  using  bubble  enhancement, 
which  allows  for  efficient  stacking  of  compound  lesions.  The 
thermal  lesion  had  a  sharp  back  wall  margin,  allowing  precise 
protection  of  tissue  structures  distal  to  the  targeted  region. 
Perfluorocarbon  emulsions  demonstrate  great  potential  for 
enhancing  clinically  relevant  HIFU  performance,  with  four- 
times  faster  treatment  times  as  well  as  improved  lesion  homo¬ 
geneity.  This  initial  phantom  study  serves  as  an  appropriate 
test  of  our  methods  and  also  provides  preliminary  indications 
for  in  vivo  animal  study  and  clinical  trials.  Efficient  lesion  for¬ 
mation  with  real-time  imaging  feedback  at  droplet  concentra¬ 
tions  less  than  1%  of  tested  safety  limits  supports  feasibility 
for  high  biocompatibility  of  this  mechanism  to  accelerate  ther¬ 
mal  therapy.  This  technology  may  be  applicable  across  a 
broad  range  of  HIFU  applications. 

ACKNOWLEDGMENTS 

The  study  was  supported  in  part  by  NIH  Grant  No. 
5R01EB000281,  DOD/BCRP  Grant  No.  BC095397P1  and 
University  of  Michigan  (UM)-Shanghai  Jiao  Tong  University 
(SJTU)  Collaboration  on  Biomedical  Technologies  grant 
award. 

Kripfgans  et  al.:  Droplet  vaporization  for  patterned  therapy  543 


Bamber,  J.  C.,  and  Hill,  C.  R.  (1979).  “Ultrasonic  attenuation  and  propaga¬ 
tion  speed  in  mammalian  tissues  as  a  function  of  temperature,” 
Ultrasound.  Med.  Biol.  5,  149-157. 

Billard,  B.  E.,  Hynynen,  K.,  and  Roemer,  R.  B.  (1990).  “Effects  of  physical 
parameters  on  high  temperature  ultrasound  hyperthermia,”  Ultrasound 
Med.  Biol.  16(4),  409^120. 

Chung,  D.  J.,  Cho,  S.  H.,  Lee,  J.  M.,  and  Hahn,  S.  T.  (2012).  “Effect  of 
microbubble  contrast  agent  during  high  intensity  focused  ultrasound  abla¬ 
tion  on  rabbit  liver  in  vivo,”  Eur.  J.  Radiol.  81(4),  e519-e523. 

Coussios,  C.,  and  Roy,  R.  (2008).  “Applications  of  acoustics  and  cavitation 
to  noninvasive  therapy  and  drug  delivery,”  Annu.  Rev.  Fluid.  Mech.  40, 
395^120. 

Crouzet,  S.,  Rebillard,  X.,  Chevallier,  D.,  Rischmann,  P.,  Pasticier,  G., 
Garcia,  G.,  Rouviere,  O.,  Chapelon,  J.  Y.,  and  Gelet,  A.  (2010). 
“Multicentric  oncologic  outcomes  of  high-intensity  focused  ultrasound  for 
localized  prostate  cancer  in  803  patients,”  Eur.  Urol.  58,  559-566. 

Fabiilli,  M.  L.,  Haworth,  K.  J.,  Nasir,  F.,  Kripfgans,  O.  D.,  Carson,  P.  L., 
and  Fowlkes,  J.  B.  (2009).  “The  role  of  inertial  cavitation  in  acoustic  drop¬ 
let  vaporization,”  IEEE  Trans.  Ultrason.  Ferroelectr.  Freq.  Control.  56, 
1006-1017. 

Fabiilli,  M.  L.,  Piert,  M.  R.,  Koeppe,  R.  A.,  Sherman,  P.  S.,  Quesada,  C.  A., 
and  Kripfgans,  O.  D.  (2013).  “Assessment  of  the  biodistribution  of  an 
[18F]FDG-loaded  perfluorocarbon  double  emulsion  using  dynamic  micro- 
PET  in  rats,”  Contrast  Media  Mol.  Imaging  8(4),  366-374. 

Frizzell,  L.  A.  (1988).  “Threshold  dosages  for  damage  to  mammalian  liver 
by  high-intensity  focused  ultrasound,”  IEEE  Trans  Ultrason  Ferroelectr 
Freq  Control.  35,  578-581. 

Gervais,  D.,  McGovern,  F.,  Arellano,  R.,  McDougal,  W.  S.,  and  Mueller, 
P.  (2003).  “Renal  cell  carcinoma:  Clinical  experience  and  technical 
success  with  radio-frequency  ablation  of  42  tumors,”  Radiology  226, 
417—424. 

Hill,  C.  R.,  and  ter  Haar,  G.  R.  (1995).  “Review  article:  High  intensity 
focused  ultrasound-potential  for  cancer  treatment,”  Br.  J.  Radiol.  68(816), 
1296-1303. 

Illing,  R.  O.,  Kennedy,  J.  E.,  Wu,  F.,  ter  Haar,  G.  R.,  Protheroe,  A.  S., 
Friend,  P.  J.,  Gleeson,  F.  V.,  Cranston,  D.  W.,  Phillips,  R.  R.,  and 
Middleton,  M.  R.  (2005).  “The  safety  and  feasibility  of  extracorporeal 
high-intensity  focused  ultrasound  (HIFU)  for  the  treatment  of  liver 
and  kidney  tumours  in  a  Western  population,”  Br.  J.  Cancer  93, 
890-895. 

Kripfgans,  O.  D.  (2002). “Acoustic  droplet  vaporization  for  diagnostic  and 
therapeutic  applications,”  Ph.D.  thesis,  University  of  Michigan,  Ann 
Arbor,  Michigan. 

Kripfgans,  O.  D.,  Fowlkes,  J.  B.,  Miller,  D.  L.,  Eldevik,  O.  P.,  and  Carson, 
P.  L.  (2000).  “Acoustic  droplet  vaporization  for  therapeutic  and  diagnostic 
applications,”  Ultrasound  Med.  Biol.  26(7),  1177-1189. 

Lambelet,  P.,  Ducret,  F.,  Leuba,  J.,  and  Geoffroy,  M.  (1991).  “Low-field  nu¬ 
clear  magnetic  resonance  relaxation  study  of  the  thermal  denaturation  of 
transferrins,”  J.  Agric.  Food  Chem.  39,  287-292. 

Li,  J.-J.,  Gu,  M.-F.,  Luo,  G.-Y.,  Liu,  L.-Z.,  Zhang,  R.,  and  Xu,  G.-L.  (2009). 
“Complications  of  high  intensity  focused  ultrasound  for  patients  with  he¬ 
patocellular  carcinoma,”  Technol.  Cancer  Res.  Treat.  8(3),  217-224. 


Lo,  A.  H.,  Kripfgans,  O.  D.,  Carson,  P.  L.,  and  Fowlkes,  J.  B.  (2006). 
“Spatial  control  of  gas  bubbles  and  their  effects  on  acoustic  fields,” 
Ultrasound  Med.  Biol.  32(1),  95-106. 

Parker,  K.  J.  (1983).  “Ultrasonic  attenuation  and  absorption  in  liver  tissue,” 
Ultrasound  Med.  Biol.  9(4),  363-369. 

Peng,  S.,  Xiong,  Y.,  Lia,  K.,  He,  M.,  Deng,  Y.,  Li  Chen,  Zou,  M.,  Chen,  W. 
Wang,  Z.,  He,  J.,  and  Zhang,  L.  (2012).  “Clinical  utility  of  a  microbubble¬ 
enhancing  contrast  (‘SonoVue’)  in  treatment  of  uterine  fibroids  with  high 
intensity  focused  ultrasound:  A  retrospective  study,”  Eur.  J.  Radiol. 
81(12),  3832-3838. 

Schwickert,  H.  C.,  Roberts,  T.  P.  L.,  Shames,  D.  M.,  van  Dijke,  C.  F., 
Disston,  A.,  Miihler,  A.,  Mann,  J.  S.,  and  Brasch,  R.  C.  (2005). 
“Quantification  of  liver  blood  volume:  Comparison  of  ultra  short  TI  inver¬ 
sion  recovery  echo  planar  imaging  (ULSTIR-EPI),  with  dynamic  30- 
gradient  recalled  echo  imaging,”  Magn  Reson  Med.  34(6),  845-852. 

Shaw,  A.,  and  Hodnett  M.  (2008).  “Calibration  and  measurement  issues  for 
therapeutic  ultrasound,”  Ultrasonics  48,  234-252. 

Stadler,  A.,  Schima,  W.,  Ba-Ssalamah,  A.,  Kettenbach,  J.,  and  Eisenhuber, 
E.  (2007).  “Artifacts  in  body  MR  imaging:  Their  appearance  and  how  to 
eliminate  them,”  Eur  Radiol.  17,  1242-1255. 

Takegami,  K.,  Kaneko,  Y.,  Watanabe,  T.,  Maruyama,  T.,  Matsumoto,  Y., 
and  Nagawa,  H.  (2004).  “Polyacrylamide  gel  containing  egg  white  as  new 
model  for  irradiation  experiments  using  focused  ultrasound,”  Ultrasound 
Med.  Biol.  30,  1419-1422. 

ter  Haar,  G.,  and  Coussios,  C.  (2007).  “High  intensity  focused  ultra¬ 
sound:  Physical  principles  and  devices,”  Int.  J.  Hyperthermia  23(2), 
89-104. 

Tung,  Y.,  Liu,  H.,  Wu,  C.,  Ju,  K.,  Chen,  W.,  and  Lin  W.  (2006).  “Contrast- 
agent-enhanced  ultrasound  thermal  ablation,”  Ultrasound  Med.  Biol.  32, 
1103-1110. 

Yamakado,  K.,  Nakatsuka,  A.,  Kobayashi,  S.,  Akeboshi,  M.,  Takaki,  H., 
Kariya,  Z.,  Kinbara,  H.,  Arima,  K.,  Yanagawa,  M.,  Hori,  Y.,  Kato,  H., 
Sugimura,  Y.,  and  Takeda,  K.  (2006).  “Radiofrequency  ablation  combined 
with  renal  arterial  embolization  for  the  treatment  of  unresectable  renal  cell 
carcinoma  larger  than  3.5  cm:  Initial  experience,”  Cardio.  Vascular  Inter. 
Radiology  29(3),  389-394. 

Zhang,  L.,  Chen,  W.-Z.,  Liu,  Y.-J.,  Hu,  X.,  Zhou,  K.,  Chen,  L.,  Peng,  S., 
Zhu,  H.,  Zou,  H.-L.,  Bai,  J.,  and  Wang,  Z.-B.  (2010b).  “Feasibility  of 
magnetic  resonance  imaging-guided  high  intensity  focused  ultrasound 
therapy  for  ablating  uterine  fibroids  in  patients  with  bowel  lies  anterior  to 
uterus,”  Eur.  J.  Radiol.  73(2),  396^103. 

Zhang,  M.,  Fabiilli,  M.  L.,  Fowlkes,  J.  B.,  Padilla,  F.,  Kripfgans,  O.  D., 
Swanson,  D.  S.,  and  Carson,  P.  L.  (2011).  “Acoustic  droplet  vaporization 
for  enhancement  of  thermal  ablation  by  high  intensity  focused  ultra¬ 
sound,”  Acad.  Rad.  18(9),  1123-1132. 

Zhang,  M.,  Fabiilli,  M.  L.,  Haworth,  K.  J.,  Fowlkes,  J.  B.,  Kripfgans,  O.  D., 
Roberts,  W.  W.,  Ives,  K.  A.,  and  Carson,  P.  L.  (2010a).  “Initial  investiga¬ 
tion  of  acoustic  droplet  vaporization  for  occlusion  in  canine  kidney,” 
Ultrasound  Med.  Biol.  36(10),  1691-1703. 

Zhang,  P.,  Kopechek,  J.  A.,  and  Porter,  T.  M.  (2013).  “The  impact  of  vapor¬ 
ized  nanoemulsions  on  ultrasound-mediated  ablation,”  J.  Therapeutic 
Ultrasound  1(2),  1-13. 


544  J.  Acoust.  Soc.  Am.,  Vol.  135,  No.  1,  January  2014 


Kripfgans  etal Droplet  vaporization  for  patterned  therapy 


IEEE  TRANSACTIONS  ON  BIOMEDICAL  ENGINEERING,  VOL.  59,  NO.  9,  SEPTEMBER  2012 


2431 


A  Preclinical  System  Prototype  for  Focused 
Microwave  Thermal  Therapy  of  the  Breast 

John  Stang*,  Member,  IEEE ,  Mark  Haynes,  Member,  IEEE ,  Paul  Carson,  Senior  Member,  IEEE , 

and  Mahta  Moghaddam,  Fellow,  IEEE 


Abstract — A  preclinical  prototype  of  a  transcutaneous  thermal 
therapy  system  has  been  developed  for  the  targeted  treatment  of 
breast  cancer  cells  using  focused  microwaves  as  an  adjuvant  to 
radiation,  chemotherapy,  and  high-intensity-focused  ultrasound. 
The  prototype  system  employs  a  2-D  array  of  tapered  microstrip 
patch  antennas  operating  at  915  MHz  to  focus  continuous-wave 
microwave  energy  transcutaneously  into  the  pendent  breast  sus¬ 
pended  in  a  coupling  medium.  Prior  imaging  studies  are  used  to 
ascertain  the  material  properties  of  the  breast  tissue,  and  these 
data  are  incorporated  into  a  multiphysics  model.  Time-reversal 
techniques  are  employed  to  find  a  solution  (relative  amplitudes  and 
phase)  for  focusing  at  a  given  location.  Modeling  tests  of  this  time- 
reversal  focusing  method  have  been  performed,  which  demonstrate 
good  targeting  accuracy  within  heterogeneous  breast  tissue.  Exper¬ 
imental  results  using  the  laboratory  prototype  to  perform  focused 
heating  in  tissue-mimicking  gelatin  phantoms  have  demonstrated 
1.5-cm-diameter  focal  spot  sizes  and  differential  heating  at  the  de¬ 
sired  focus  sufficient  to  achieve  an  antitumor  effect  confined  to  the 
target  region. 

Index  Terms — Hyperthermia,  microstrip  patch  antennas,  mi¬ 
crowave  imaging,  microwave  therapy,  time  reversal. 

I.  Introduction 

THERMAL  therapies — particularly  the  use  of  thermal  ab¬ 
lation  (including  microwave,  RF,  cryoablation,  and  high- 
intensity-focused  ultrasound  [1]),  hyperthermia,  and  heat  acti¬ 
vated  drug  delivery  in  the  treatment  of  cancer — have  been  the 
focus  of  increasing  laboratory  and  clinical  research.  In  the  vari¬ 
ous  ablation  methods,  either  significantly  elevated  temperatures 
(typically  above  50  °C)  generated  by  electromagnetic  waves  or 
ultrasound,  or  freezing  temperatures  generated  by  cryogens  are 
used  to  cause  very  rapid  and  localized  tissue  destruction.  In  the 
case  of  hyperthermia,  moderately  elevated  temperatures  (typi¬ 
cally  between  40  °C  and  45  °C)  have  been  shown  to  achieve 
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Fig.  1.  Focused  microwave  therapy  of  the  breast:  General  approach. 


cytotoxic  effects  that  render  cancer  cells  more  vulnerable  to 
radiotherapy  [2]  and  chemotherapy  [3],  as  well  as  inducing 
both  apoptotic  and  necrotic  cell  death  given  a  sufficient  thermal 
dose  [4]-[6]. 

Motivated  by  this  knowledge,  researchers  have  developed  mi¬ 
crowave  systems  capable  of  inducing  localized  hyperthermia  as 
an  adjuvant  cancer  treatment  modality.  Laboratory  and  clini¬ 
cal  studies  of  a  system  using  a  two-element  adaptive  phased 
array  and  an  invasive  electric-field  probe  to  focus  within  a  com¬ 
pressed  breast  have  been  reported  by  Dooley  et  al.  [7]  and 
Vargas  et  al  [8].  In  addition,  Stauffer  et  al  [9],  [10]  and  De- 
whirst  et  al  [11]  have  demonstrated  MR-guided  phased  arrays 
for  the  treatment  of  locally  advanced  breast  cancer,  chest  wall 
recurrence,  and  osteosarcomas.  Several  groups  have  also  be¬ 
gun  clinical  testing  of  temperature-triggered  drug  delivery  using 
low  temperature- sensitive  liposomes  activated  with  hyperther¬ 
mia  [12]  and  RF  ablation  [13].  A  summary  of  the  methods  and 
recent  progress  in  heat-activated  drug  delivery  can  be  found 
in  [14]. 

With  these  potential  applications  in  mind,  a  preclinical  pro¬ 
totype  of  a  noninvasive  thermal  therapy  system  has  been  de¬ 
veloped  for  the  targeted  treatment  of  breast  cancer  cells  using 
focused  microwaves.  In  this  system,  an  array  of  antennas  oper¬ 
ating  at  915  MHz  is  used  to  focus  continuous- wave  microwave 
energy  transcutaneously  into  the  pendent  breast  suspended  in 
a  coupling  medium  (shown  in  Fig.  1).  Prior  imaging  studies 
are  used  to  ascertain  the  material  properties  of  the  breast  tis¬ 
sue,  and  these  data  are  incorporated  into  a  multiphysics  model. 
Time-reversal  techniques  (described  in  Section  II-A)  are  then 
employed  to  find  a  solution  (relative  amplitudes  and  phase)  for 
focusing  at  a  given  location,  resulting  in  maximal  thermal  dose 
at  the  tumor  location.  Using  this  system,  focal  spot  sizes  in 
the  array  plane  of  approximately  1.5  cm  in  diameter  have  been 
achieved  (discussed  in  Section  III-B),  and  significant  differen¬ 
tial  heating  in  the  target  region  has  been  observed  in  focused 
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heating  tests  within  tissue-mimicking  gelatin  phantoms  (details 
in  Section  III-C). 

Based  on  these  results,  this  system  has  the  potential  to  offer 
improved  targeting  and  delivery  of  focused  heating  over  current 
microwave  thermal  therapy  systems,  without  the  need  for  inva¬ 
sive  probes.  Upon  completion  of  an  optimized  clinical  version, 
this  system  will  have  the  potential  to  be  used  in  neoadjuvant 
thermal  treatment  for  presurgical  tumor  reduction,  noninvasive 
tumor  ablation,  locally  enhanced  drug  delivery,  and  in  postop¬ 
erative  adjuvant  therapy  to  reduce  local  recurrence. 

II.  Methods 

Development  of  the  preclinical  microwave  thermal  therapy 
prototype  system  was  completed  in  a  series  of  stages.  First,  an 
initial  2-D  antenna  array  was  designed,  and  a  full-wave  forward 
model  of  the  microwave  therapy  array  was  created.  Following 
that,  modeling  tests  of  the  time-reversal  focusing  method  were 
performed  to  ensure  good  targeting  accuracy  within  heteroge¬ 
neous  breast  tissue.  Next,  the  focusing  ability  of  the  system 
was  characterized  by  measuring  the  microwave  power  distribu¬ 
tion  within  the  therapy  chamber  using  a  coaxial  probe.  Finally, 
tissue-mimicking  gelatin  phantoms  were  developed  and  exper¬ 
iments  were  performed  to  test  differential  heating  capability  at 
the  desired  focus,  along  with  the  system’s  ability  to  achieve  a 
tumoricidal  thermal  dose  confined  to  the  target  region. 

A.  Image-Based  Time-Reversal  Focusing 

An  early  study  by  Surowiec  and  Bicher  [15]  proposed  a  three- 
element-focused  hyperthermia  system  that  identified  the  need 
to  account  for  the  various  effects  of  source  and  array  geometry 
on  focusing.  Since  that  study,  researchers  have  developed  var¬ 
ious  methods  of  computing  the  relative  magnitude  and  phase 
parameters  necessary  to  focus  microwave  energy  at  the  target. 
In  the  clinical  system  reported  by  the  authors  in  [16]  and  [17], 
an  invasive  electric-field  probe  placed  within  the  compressed 
breast  is  used  to  focus  a  two-element  adaptive  phased  array.  Re¬ 
cent  theoretical  studies  have  also  reported  less  invasive  methods 
including  the  use  of  a  deformable  mirror  method  [18]  and  trans¬ 
mit  beamforming  methods  [19].  Finally,  there  are  time-reversal 
techniques  which  are  rooted  in  theory  previously  developed 
for  ultrasound  focusing  [20],  such  as  those  discussed  in  [21] 
and  [22]. 

In  the  patient- specific  variant  of  the  time-reversal  method 
employed  here,  the  predicted  microwave  tissue  properties  from 
a  prior  imaging  study  (e.g.,  ultrasound,  MRI,  CT,  or  possibly 
microwave  imaging)  are  incorporated  into  a  full- wave  forward 
model  of  the  microwave  therapy  array  [see  Fig.  2(a)] .  A  synthetic 
dipole  point  source  polarized  along  the  axis  of  the  cylindrical 
cavity  is  then  excited  at  the  target  location  within  the  breast  [as 
shown  in  the  field  distribution  on  the  left  in  Fig.  2(b)].  The  rela¬ 
tive  phase  of  the  electric  field  at  each  antenna  feed  location  after 
traveling  from  the  point  source  excitation  through  the  interven¬ 
ing  breast  tissue  and  coupling  fluid  is  computed  and  stored.  The 
complex  phase  information  is  then  conjugated  (time-reversed) 
to  produce  the  phase  delays  needed  to  electronically  steer  the 
focus  of  the  microwave  antenna  array  back  to  the  inclusion  loca¬ 
tion  [as  shown  in  the  field  distribution  on  the  right  in  Fig.  2(b)]. 


(a)  (b) 


Fig.  2.  (a)  Sample  imported  MR  image  of  fatty  breast,  (b)  Patient- specific 

time-reversal  focusing  using  imported  MR  image.  (Left)  Outgoing  field  distri¬ 
bution  due  to  synthetic  dipole  placed  at  the  center  of  the  array.  (Right)  Incoming 
field  distribution  focused  using  time-reversed  phase  data. 


Fig.  3.  Two-dimensional  microwave  thermal  therapy  system  prototype, 
(a)  Signal  generator,  (b)  12-channel  electronic  phase  control,  (c)  10-W  per 
channel  amplifiers  (120  W  total),  and  (d)  12-element  therapy  array. 

A  table  of  phase  delays  needed  to  focus  within  the  therapy  cham¬ 
ber  can  then  be  generated  by  repeating  this  procedure  for  the  set 
of  desired  focal  locations.  In  the  experiments  that  follow,  this 
was  done  on  a  1.5-cm  uniform  grid  of  49  possible  focal  states. 

B.  System  Design 

As  shown  in  Fig.  3,  the  prototype  microwave  therapy  system 
consists  of  a  signal  generator,  a  12-channel  power  divider  and 
electronically  controlled  phase  shifter  network,  12  microwave 
amplifiers  providing  10  W  of  power  per  channel,  and  a  15-cm- 
diameter  cylindrical  cavity  containing  a  12-element  antenna 
array.  The  signal  generator  provides  a  low-power  915-MHz 
continuous  microwave  source,  which  is  first  equally  divided 
to  each  channel  using  the  Minicircuits  power  divider  (ZN12PD- 
252-S+).  Each  channel  is  then  simultaneously  phase-delayed 
using  a  series  of  two  electronically  controlled  phase  shifters 
(Minicircuits  JSPH-1000)  by  an  amount  determined  with  the 
time-reversal  focusing  method  discussed  previously.  Each  of  the 
appropriately  delayed  signals  are  then  amplified  with  a  10-W 
Minicircuits  broadband  amplifier  (ZHL-10W-2G+)  and  trans¬ 
mitted  into  the  therapy  chamber  using  a  microstrip  patch  antenna 
designed  to  operate  with  maximum  efficiency  at  915  MHz.  The 
loss  through  the  power  divider/phase  shifter  network  was  mea¬ 
sured  to  be  14.5  dB  and  is  compensated  for  using  the  signal 
generator.  The  power  at  the  output  of  each  amplifier  was  mea¬ 
sured  to  be  40.7  dB-m  (1 1.75  W)  using  a  spectrum  analyzer. 
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(a)  (b) 


Fig.  4.  (a)  Initial  temperature  measurements  with  digital  probe  thermome¬ 

ters  in  50/50  oil-in-water  gelatin  phantom  with  30/70  oil-in-water  cylindrical 
inclusion,  (b)  Thermistor  array  above  homogeneous  50/50  oil-in-water  gelatin 
phantom  in  60/40  matching  fluid. 


Fig.  5.  Coaxial  probe  for  mapping  microwave  power  distribution  in  the  plane 
of  the  therapy  array  (shown  in  photo  without  coupling  fluid  for  clarity). 


1 )  Microstrip  Antenna  Array:  In  microwave  thermal  ther¬ 
apy,  as  in  near-field  microwave  imaging,  antennas  must  be  de¬ 
signed  to  radiate  either  directly  into  tissue,  or  into  a  matching 
fluid  medium  designed  to  mimic  the  dielectric  properties  of  the 
tissue.  Also,  mutual  impedance  effects  can  negatively  impact 
antenna  performance  in  a  variety  of  ways  (shift  resonant  fre¬ 
quency,  alter  radiation  pattern,  reduce  matching  performance). 
Furthermore,  these  effects  can  influence  each  antenna  in  an  ar¬ 
ray  differently,  depending  on  the  geometry  of  the  antenna  and 
the  array.  Since  it  is  not  feasible  to  account  for  all  of  these  effects 
analytically,  the  general  strategy  was  to  use  an  iterative  combi¬ 
nation  of  analytical  methods  and  computer  simulation  in  order 
to  design  the  proof  of  concept  array.  First,  a  standard  model  of 
a  rectangular  patch  antenna  was  used  to  make  an  initial  design 
intended  to  resonate  at  915  MHz  (ISM  Band).  Then,  the  geom¬ 
etry  was  modified  with  a  linear  taper  to  account  for  the  effect 
of  the  dielectric  loading  of  the  biological  tissue  on  the  input 
impedance.  This  taper  was  subsequently  optimized  through  an 
iterative  process  of  simulations  guided  by  the  analytical  model. 
Finally,  the  array  geometry  and  antenna  spacing  were  optimized 
using  a  similar  iterative  process  of  computer  simulation.  Further 
details  can  be  found  in  [23]  and  [24]. 

2)  Emulsion  Matching  Medium  and  Gelatin  Phantoms:  A 
surfactant-stabilized  two-phase  oil-in-water  emulsion  is  used 
as  a  matching  medium  designed  for  optimal  coupling  of  mi¬ 
crowave  power  between  the  array  and  the  breast.  Oil-in- water 
emulsions  were  chosen  due  to  their  straightforward  fabrication 
from  readily  available,  inexpensive  materials,  their  relatively 
low  attenuation,  low  viscosity  for  ease  of  circulation,  and  the 
ability  to  adjust  their  dielectric  properties  to  match  the  tissue  of 
interest  simply  by  varying  the  ratio  of  oil  and  water.  The  emul¬ 
sion  coupling  medium  used  in  the  current  therapy  system  is 
composed  of  60%  vegetable  oil,  36%  DI  water,  and  4%  HLB10 
surfactant  (46%  SPAN  80  and  54%  TWEEN  80)  and  was  emul¬ 
sified  using  an  ultrasonic  sonicator.  The  dielectric  properties  of 
the  emulsion  coupling  medium  were  measured  using  the  Agilent 
85070E  open-ended  coaxial  dielectric  probe  kit  to  be  er  =22.9 
and  a  =  0.07  at  915  MHz.  In  order  to  dissipate  the  heat  gener¬ 
ated  by  the  absorption  of  the  significant  electromagnetic  energy 
present  in  the  immediate  near  field  of  the  antenna  elements, 
the  emulsion  matching  medium  was  circulated  and  cooled  us¬ 
ing  an  ice  water  bath.  Both  heterogeneous  [see  Fig.  4(a)]  and 
homogeneous  [see  Fig.  4(b)]  solid  phantoms  were  created  with 
a  range  of  dielectric  properties  by  adding  unflavored  gelatin  to 


emulsions  with  varying  ratios  of  oil/water  that  were  heated  to 
slightly  below  boiling.  Further  details  of  the  emulsion  coupling 
fluid  and  tissue-mimicking  phantoms  can  be  found  in  [25]. 

C.  System  Characterization 

1 )  Characterization  of  Microwave  Power  Delivery:  In  order 
to  evaluate  the  focusing  ability  of  the  system,  the  microwave 
power  distribution  within  the  therapy  chamber  filled  with  ho¬ 
mogeneous  coupling  fluid  was  measured  using  a  coaxial  probe 
(shown  in  Fig.  5).  The  probe  was  manually  positioned  on  a 
2-D  grid  with  5 -mm  spacing,  and  the  field  for  49  electronically 
scanned  focal  states  was  measured  at  each  position.  From  these 
data,  microwave  power  distribution  maps  were  created  for  each 
focal  state,  and  the  results  are  discussed  in  Section  III-B. 

2 )  Characterization  of  Thermal  Distribution:  A  series  of  fo¬ 
cused  heating  experiments  was  conducted  on  a  number  of  tissue- 
mimicking  gelatin  phantoms.  To  initially  confirm  that  focused 
heating  was  occurring  at  the  desired  location,  temperature  mea¬ 
surements  were  taken  after  a  series  of  heating  times  using  simple 
digital  probe  thermometers  at  several  sample  locations  within 
a  gelatin  phantom  and  surrounding  emulsion  matching  medium 
[shown  in  Fig.  4(a)].  Following  these  simple  temperature  probe 
experiments,  a  linear  array  of  thermistors  [as  shown  in  Fig.  4(b)] 
was  used  to  map  the  temperature  distribution  throughout  the 
therapy  array  for  a  number  of  gelatin  phantoms,  focal  locations, 
and  heating  times.  This  was  accomplished  by  first  orienting  the 
seven-element  thermistor  array  along  one  of  the  principal  axes 
of  the  cavity.  The  array  was  continually  lowered  by  hand  to  20 
points  in  depth  at  0.5 -cm  increments,  where  the  temperature 
measured  by  each  thermistor  was  recorded.  This  was  done  for 
seven  different  positions  across  the  cavity,  spaced  at  1.5-cm  in¬ 
crements,  in  order  to  create  a  7  x  7  x  20  array  of  temperature 
samples  which  span  the  volume  of  the  gelatin  phantoms  as  well 
as  a  portion  of  the  surrounding  coupling  medium. 

While  temperature  probes  are  used  here,  in  the  clinical  setting, 
an  integrated  noninvasive  method  of  real-time  thermal  monitor¬ 
ing  will  be  needed.  One  possible  method  would  involve  using  a 
ring  array  of  transducer  elements  integrated  into  the  current  mi¬ 
crowave  array  to  perform  ultrasonic  transmission  tomography. 
Using  this  method,  the  opposite  sign  of  the  thermal  dependence 
of  speed  of  sound  in  fat-  and  water-bearing  tissues  can  be  ac¬ 
counted  for  given  the  substantially  different  speeds  of  sound  of 
these  two  tissue  classes  [26],  [27].  Imaging  temperature  with 
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Return  Loss 


Fig.  6.  HFSS  predicted  versus  measured  antenna  reflection  coefficients. 


TABLE  I 

Dielectric  Values  Used  in  Numerical  Breast  Phantoms 


Material 

er 

a  (S/m) 

Glandular  Tissue  (1.1,  1.2,  1.3) 
Transitional  Tissue  (2) 
Adipose  Tissue  (3.1,  3.2,  3.3) 
Skin 

Muscle  (Chest  Wall) 
Matching  Fluid 

55,  45,  35 
25 

15,  10,  5 

40 

65 

24 

0.6,  0.5,  0.4 
0.3 

0.2,  0.15,  0.1 
1.0 

2.0 

0.5 

smaller  aperture  ultrasound  is  also  possible,  e.g.,  [28]— [30],  as 
is  microwave  radiometry  [31].  While  real-time  absolute  temper¬ 
ature  measurement  is  desirable,  ultrasonic  or  other  imaging  of  a 
final  endpoint,  such  as  elasticity  changes  from  thermal  necrosis, 
may  be  sufficient  as  well. 

III.  Results  and  Discussion 

A.  Antenna  Performance 

The  tapered  microstrip  antennas  used  in  the  therapy  array 
were  designed  using  Ansoft  HFSS  to  have  optimal  efficiency 
when  radiating  at  915  MHz  into  a  coupling  medium  with  di¬ 
electric  properties  of  er  =  23  and  a  =  0.06.  A  comparison  of 
the  HFSS  predicted  and  measured  reflection  coefficients  for  each 
of  the  12  therapy  antennas  is  shown  in  Fig.  6,  demonstrating  that 
all  of  the  antennas  are  well  matched  (reflection  coefficient  less 
than  — 15  dB)  at  915  MHz  radiating  into  the  coupling  medium. 

B.  Array  Focusing 

1 )  Synthetic  Focusing  Tests  in  Anatomically  Realistic  Numer¬ 
ical  Breast  Phantoms:  The  focusing  capabilities  of  the  time- 
reversal  method  were  first  tested  synthetically  using  a  ten- 
element  array  to  focus  within  a  variety  of  anatomically  real¬ 
istic  numerical  breast  phantoms.  These  modeling  tests  were 
performed  by  importing  MRI-based  numerical  breast  phan¬ 
toms  [32]  that  spanned  the  range  of  BIRADS  breast  types  [33] 
into  COMSOL,  which  was  then  used  to  simulate  the  interac¬ 
tion  of  the  time-reversal-focused  electromagnetic  waves  with 
the  reported  [34],  [35]  electrical  properties  (see  Table  I)  of 
the  various  breast  tissue  types.  These  results  (shown  in  Fig.  7) 


demonstrate  the  ability  of  the  system  to  achieve  focusing  within 
a  1.5 -cm  target  region — for  breasts  of  widely  differing  density 
and  heterogeneity — without  prohibitive  microwave  power  de¬ 
livery  to  other  regions  of  the  breast. 

2 )  Experimental  Focusing  Tests  in  Coupling  Media:  The  mi¬ 
crowave  power  distribution  for  49  focal  states  was  measured  us¬ 
ing  a  coaxial  probe  that  was  manually  positioned  on  a  2-D  grid 
with  5 -mm  spacing  (see  Fig.  5).  Microwave  power  distribution 
maps  for  nine  sample  focal  states  are  shown  in  Fig.  8.  These 
results  suggest  that  using  the  current  time-reversal  method  with 
a  12-element,  15-cm-diameter  cylindrical  array  (each  element 
7.5  cm  from  the  center)  can  achieve  suitable  focusing  within  a 
circular  region  approximately  10  cm  in  diameter.  This  is  evi¬ 
dent  by  good  targeting  accuracy,  focal  spot  size  and  shape,  and 
power  delivered  to  the  target  location  in  this  region.  However, 
near  the  perimeter  of  the  chamber,  we  observe  reduced  targeting 
accuracy,  deformed  and  reduced  focal  spots,  and  significant  side 
lobes,  which  are  expected.  These  limitations  can  be  overcome 
with  additional  elements  or  modifications  to  the  geometry  of  the 
array  as  needed  for  a  given  focusing  requirement. 

C.  Thermal  Distribution  in  Tissue-Mimicking  Phantoms 

In  the  first  series  of  phantom  heating  experiments,  tempera¬ 
ture  measurements  were  taken  using  digital  probe  thermometers 
at  sample  locations  within  a  gelatin  phantom  and  surrounding 
emulsion  matching  medium  [see  Fig.  4(a)].  Following  these  sim¬ 
ple  temperature  probe  experiments,  a  linear  array  of  thermistors 
[see  Fig.  4(b)]  was  used  to  map  the  temperature  distribution 
throughout  the  phantom  and  surrounding  fluid  for  a  number  of 
gelatin  phantoms,  focal  locations,  and  heating  times.  The  re¬ 
sults  of  both  the  thermistor  array  mapping  and  the  digital  probe 
measurements  for  a  focal  spot  offset  3  cm  from  the  center  of 
the  therapy  chamber  are  presented  in  Fig.  9.  In  these  tests,  as  in 
others  performed  on  phantoms  of  different  sizes  and  composi¬ 
tion  (some  containing  cancer  mimicking  inclusions),  significant 
differential  heating  was  observed  at  the  target  location. 

Finally,  temperature  maps  of  the  focal  planes  are  shown  for 
both  a  3-cm  offset  case  (see  Fig.  10)  and  a  center-focused 
case  (see  Fig.  11).  These  results  show  that  the  temperature 
distribution  is  in  good  agreement  with  the  earlier  microwave 
power  maps.  The  asymmetry  around  the  focal  spots  along  the 
x-direction  seen  in  the  offset  focus  case  is  most  likely  due  to 
imprecision  in  the  manual  placement  of  the  thermistor  array. 
Despite  their  somewhat  low  resolution,  they  confirm  the  ability 
of  the  system  to  achieve  an  antitumor  effect  within  a  1.5 -cm 
target  region  in  the  array  plane  without  unwanted  heating  in  any 
other  region. 

IV.  Conclusion 

A  preclinical  prototype  of  a  microwave  thermal  therapy  sys¬ 
tem  that  employs  an  image-based  time-reversal  method  has  been 
developed  and  tested.  Simulated  and  experimental  results  have 
demonstrated  that  such  a  system  can  achieve  suitable  focusing 
and  heating  for  the  potential  application  of  thermal  therapy  of 
the  breast.  Based  on  the  results  of  these  studies,  the  design  and 
fabrication  of  a  clinically  viable  focused  microwave  array  is  cur¬ 
rently  underway  that  will  demonstrate  full  3-D  focusing  ability 
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Fatty  Breast:  Dielectric  Map 


Coronal 


Transverse 


Time-Reversal  Focusing 


Coronal  Transverse 


Fig.  7.  COMSOL  predicted  microwave  power  distributions  for  time-reversal  focusing  in  (a)  fatty  and  (b)  dense  breast  (ten-element  array).  Elongation  in  transverse 

plane  due  to  2-D  array  focusing. 
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(c) 


Fig.  8.  Microwave  power  distributions  on  a  subgrid  with  1.5-cm  spacing: 
simulated  results  on  top  versus  measured  results  on  bottom,  (a)  x  =  —4.5  cm, 
(b)  x  =  —3.0  cm,  and  (c)  x  =  —1.5  cm. 

using  nonlinear  array  optimization  techniques  with  substantially 
improved  anterior  to  posterior  concentration  of  heating  and  more 
efficient  cooling  of  the  coupling  fluid.  Upon  completion,  the  fo¬ 
cused  microwave  thermal  therapy  system  will  have  the  potential 
to  be  used  in  neoadjuvant  tumor  treatment  prior  to  surgery,  to 


(b) 
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Offset  Focus  in  Homogeneous  Phantom,  Long  Heating  Times 


Fig.  9.  Focus  at  x  =  3  cm  and  y  =  0  cm.  (a)  Simulated  microwave  power 
distribution,  (b)  Measured  microwave  power  distribution,  (c)  Temperature  after 
45  min  of  heating,  measured  with  thermistor  array,  (d)  Temperature  as  a  function 
of  time,  measured  with  digital  probe  thermometers. 
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Fig.  10.  Temperature  maps  of  focused  heating  at  x  =  3  cm  and  y  =  0  cm  in  homogeneous  phantom  for  30-  and  45-min  heating  times.  Sagittal  and  transverse 
slices  show  vertical  extent  of  the  active  antennas;  coronal  slice  outlines  the  chamber. 
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Fig.  11.  Temperature  maps  of  focused  heating  at  x  =  0  cm  and  y  =  0  cm  in  homogeneous  phantom  for  30-  and  45-min  heating  times.  Sagittal  and  transverse 
slices  show  vertical  extent  of  the  active  antennas;  coronal  slice  outlines  the  chamber. 
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completely  ablate  the  tumor,  to  selectively  enhance  drug  deliv¬ 
ery  and  absorption  through  focused  thermal  release,  or  as  an 
adjuvant  to  postoperative  radiation  therapy  and  chemotherapy. 
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Malignant  tissues  have  electrical  and  acoustical  properties  that  make  them  distinguishable  from 
healthy  tissue,  including  an  ability  to  highly  absorb  microwave  and  ultrasound  energy  and  thus 
be  heated  more  quickly  than  surrounding  tissue.  This  property  can  be  exploited  for  thermal 
therapies  in  which  elevated  temperatures  are  used  to  achieve  cell  death  or  render  the  cells  more 
vulnerable  to  ionizing  radiation  and  chemotherapy.  Both  microwave  and  ultrasound  modalities 
have  been  used  for  thermal  therapies;  however,  there  are  limitations  with  each  method.  Drs 
Moghaddam  and  Carson  arc  investigating  the  combination  of  two  thermal  therapies:  high 
p>edkted  intensity  focused  ultrasound  (H1FU),  which  is  limited  in  its 

efficiency  to  treat  large  tumors  fast  with  uniformity,  and  high 
intensity  focused  microwave  (H1FW),  which  can  treat  fast 
and  relatively  uniformly  but  without  the  resolution  to  avoid 
heating  adjacent  healthy  tissues  or  to  achieve  uniform  heating 
at  sharp  boundaries  of  tissue 
properties.  The  proposed 
combined  system  seeks  to 
remove  the  limitations 
inherent  in  individual 
application  of  these 
therapies. 

Fig  la.  Predicted  field 
intensities  heside  measured 
thermal  maps  in  the  20cm 

'  'l  •'  f  ^  i  ?  b  $  W  treatment  chamber.  The 
steering  ability  of  the  array  is  clearly  demonstrated. 

Dr.  Moghaddam  is  working  to  improve  targeting  the 
microwaves  in  HJFW  to  overcome  focus  and  uniformity  issues 
using  an  array  design  that  could  be  fine-tuned  for  each  patient. 

(Fig  la).  Dr.  Carson  is  developing  a  transducer  array  that  will 
allow  placement  of  the  HJFU  focal  zone  close  to  the  chest  wall 
and  permit  ablation  of  deep  tumors  at  higher  speed  and 
uniformity  than  currently  possible  with  ultrasound  (Fig.  lb). 

Using  information  gathered  from  the  study  of  each  treatment 
modality,  the  collaborators  aim  to  combine  these  therapies  in  a 
microwave-ultrasound  synergistic  thermal  (MUST)  treatment 
system.  Experimental  results  using  laboratory  prototypes 
generated  focused  heating  sufficient  for  an  antitumor  effect  in 
tissue-mimicking  gelatin  phantoms.  The  heating  focal  spot 
sizes  were  as  small  as  1.5  cm  with  microwaves  and  1.5  mm  with 
ultrasound.  Further  development  of  the  MUST  system  will 
capitalize  on  the  strengths  of  each  modality  such  that  the 
combined  system  ould  substantially  enhance  treatment  results 
through  precise  and  uniform  heating  accompanied  by 
simultaneous  imaging  of  the  temperature. 
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Fig.  lb.  T 2 -weighted  MR 
images  of  ultrasound  thermal 
lesions  in  a  uniform  gel. 
placed  in  a  spiral  without 
(top)  and  with  (bottom) 
controlled  gas  huhblc 
gene  ra  t i on  from  1 1 r 
adm  inistered  perfluor- 
ocarhon  microdroplets .  The 
latter  greatly  enhances 
heating  only  in  the  focal  zone. 
Note  the  small,  sparse  lesions 
in  the  top,  “unenhanc  ed  " 
image  and  the  uniform 
coverage  from  larger  lesions 
in  the  “ enhanced  "  image. 
Both  were  obtained  in  the 
same  exposure  time  and 
ultrasound  output  levels . 


